# Source code for pointpats.centrography

    """
    Centrographic measures for point patterns

    TODO

    - testing
    - documentation

    """

    __author__ = "Serge Rey sjsrey@gmail.com"

    __all__ = [
        "mbr",
        "hull",
        "mean_center",
        "weighted_mean_center",
        "manhattan_median",
        "std_distance",
        "euclidean_median",
        "ellipse",
        "minimum_rotated_rectangle",
        "minimum_bounding_rectangle",
        "skyum",
        "dtot",
        "_circle",
    ]

    import sys
    import numpy as np
    import warnings
    import copy
    import math

    from math import pi as PI
    from scipy.spatial import ConvexHull
    from libpysal.cg import get_angle_between, Ray, is_clockwise
    from scipy.spatial import distance as dist
    from scipy.optimize import minimize
    import shapely

    not_clockwise = lambda x: not is_clockwise(x)

    MAXD = sys.float_info.max
    MIND = sys.float_info.min

    def minimum_bounding_rectangle(points):
        """
        Find minimum bounding rectangle of a point array.

        Parameters
        ----------
        points : arraylike
                 (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        min_x  : float
                 leftmost value of the vertices of minimum bounding rectangle.
        min_y  : float
                 downmost value of the vertices of minimum bounding rectangle.
        max_x  : float
                 rightmost value of the vertices of minimum bounding rectangle.
        max_y  : float
                 upmost value of the vertices of minimum bounding rectangle.

        """
        points = np.asarray(points)
        min_x = min_y = MAXD
        max_x = max_y = MIND
        x, y = zip(*points)
        min_x = min(x)
        min_y = min(y)
        max_x = max(x)
        max_y = max(y)
        return min_x, min_y, max_x, max_y

    def minimum_rotated_rectangle(points, return_angle=False):
        """
        Compute the minimum rotated rectangle for an input point set.

        This is the smallest enclosing rectangle (possibly rotated)
        for the input point set. It is computed using Shapely.

        Parameters
        ----------
        points : numpy.ndarray
            A numpy array of shape (n_observations, 2) containing the point
            locations to compute the rotated rectangle
        return_angle : bool
            whether to return the angle (in degrees) of the angle between
            the horizontal axis of the rectanle and the first side (i.e. length).

        Returns
        -------
        an numpy.ndarray of shape (4, 2) containing the coordinates
        of the minimum rotated rectangle. If return_angle is True,
        also return the angle (in degrees) of the rotated rectangle.

        """
        points = np.asarray(points)
        out_points = shapely.get_coordinates(
            shapely.minimum_rotated_rectangle(shapely.multipoints(points))
        )[:-1]
        if return_angle:
            angle = (
                math.degrees(
                    math.atan2(
                        out_points[1][1] - out_points[0][1],
                        out_points[1][0] - out_points[0][0],
                    )
                )
                % 90
            )
            return (out_points[::-1], angle)
        return out_points[::-1]

    def mbr(points):
        warnings.warn(
            "This function will be deprecated in the next release of pointpats.",
            FutureWarning,
            stacklevel=2,
        )
        return minimum_bounding_rectangle(points)

    mbr.__doc__ = minimum_bounding_rectangle.__doc__

    def hull(points):
        """
        Find convex hull of a point array.

        Parameters
        ----------
        points: arraylike
                (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        _     : array
                (h,2), points defining the hull in counterclockwise order.
        """

        points = np.asarray(points)
        h = ConvexHull(points)
        return points[h.vertices]

    def mean_center(points):
        """
        Find mean center of a point array.

        Parameters
        ----------
        points: arraylike
                (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        _     : array
                (2,), (x,y) coordinates of the mean center.
        """

        points = np.asarray(points)
        return points.mean(axis=0)

    def weighted_mean_center(points, weights):
        """
        Find weighted mean center of a marked point pattern.

        Parameters
        ----------
        points  : arraylike
                  (n,2), (x,y) coordinates of a series of event points.
        weights : arraylike
                  a series of attribute values of length n.

        Returns
        -------
        _      : array
                 (2,), (x,y) coordinates of the weighted mean center.
        """

        points, weights = np.asarray(points), np.asarray(weights)
        w = weights * 1.0 / weights.sum()
        w.shape = (1, len(points))
        return np.dot(w, points)[0]

    def manhattan_median(points):
        """
        Find manhattan median of a point array.

        Parameters
        ----------
        points  : arraylike
                  (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        _      : array
                 (2,), (x,y) coordinates of the manhattan median.
        """

        points = np.asarray(points)
        if not len(points) % 2:
            s = "Manhattan Median is not unique for even point patterns."
            warnings.warn(s)
        return np.median(points, axis=0)

    def std_distance(points):
        """
        Calculate standard distance of a point array.

        Parameters
        ----------
        points  : arraylike
                  (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        _      : float
                 standard distance.
        """

        points = np.asarray(points)
        n, p = points.shape
        m = points.mean(axis=0)
        return np.sqrt(((points * points).sum(axis=0) / n - m * m).sum())

    def ellipse(points):
        """
        Calculate parameters of standard deviational ellipse for a point pattern.

        Parameters
        ----------
        points : arraylike
                 (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        _      : float
                 semi-major axis.
        _      : float
                 semi-minor axis.
        theta  : float
                 clockwise rotation angle of the ellipse.

        Notes
        -----
        Implements approach from:

        https://www.icpsr.umich.edu/CrimeStat/files/CrimeStatChapter.4.pdf
        """

        points = np.asarray(points)
        n, k = points.shape
        x = points[:, 0]
        y = points[:, 1]
        xd = x - x.mean()
        yd = y - y.mean()
        xss = (xd * xd).sum()
        yss = (yd * yd).sum()
        cv = (xd * yd).sum()
        num = (xss - yss) + np.sqrt((xss - yss) ** 2 + 4 * (cv) ** 2)
        den = 2 * cv
        theta = np.arctan(num / den)
        cos_theta = np.cos(theta)
        sin_theta = np.sin(theta)
        n_2 = n - 2
        sd_x = (2 * (xd * cos_theta - yd * sin_theta) ** 2).sum() / n_2
        sd_y = (2 * (xd * sin_theta - yd * cos_theta) ** 2).sum() / n_2
        return np.sqrt(sd_x), np.sqrt(sd_y), theta

    def dtot(coord, points):
        """
        Sum of Euclidean distances between event points and a selected point.

        Parameters
        ----------
        coord   : arraylike
                  (x,y) coordinates of a point.
        points  : arraylike
                  (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        d       : float
                  sum of Euclidean distances.

        """
        points = np.asarray(points)
        xd = points[:, 0] - coord[0]
        yd = points[:, 1] - coord[1]
        d = np.sqrt(xd * xd + yd * yd).sum()
        return d

    def euclidean_median(points):
        """
        Calculate the Euclidean median for a point pattern.

        Parameters
        ----------
        points: arraylike
                (n,2), (x,y) coordinates of a series of event points.

        Returns
        -------
        _     : array
                (2,), (x,y) coordinates of the Euclidean median.

        """
        points = np.asarray(points)
        start = mean_center(points)
        res = minimize(dtot, start, args=(points,))
        return res["x"]

    def minimum_bounding_circle(points):
        """
        Implements Skyum (1990)'s algorithm for the minimum bounding circle in R^2.

        Store points clockwise.
        Find p in S that maximizes angle(prec(p), p, succ(p) THEN radius(prec(
        p), p, succ(p)). This is also called the lexicographic maximum, and is the last
        entry of a list of (radius, angle) in lexicographical order.

        * If angle(prec(p), p, succ(p)) <= 90 degrees, then finish.

        * If not, remove p from set.

        Parameters
        ----------
        points  :   numpy.ndarray
            a numpy array of shape (n_observations, 2) to compute
            the minimum bounding circle

        Returns
        -------
        (x,y),center for the minimum bounding circle.
        """
        points = hull(points)
        if not_clockwise(points):
            points = points[::-1]
            if not_clockwise(points):
                raise Exception("Points are neither clockwise nor counterclockwise")
        POINTS = copy.deepcopy(points)
        removed = []
        i = 0
        if HAS_NUMBA:
            circ = _skyum_numba(POINTS)[0]
        else:
            circ = _skyum_lists(POINTS)[0]
        return (circ[1], circ[2]), circ[0]

    def skyum(points):
        warnings.warn(
            "This function will be deprecated in the next release of pointpats.",
            FutureWarning,
            stacklevel=2,
        )
        return minimum_bounding_circle(points)

    skyum.__doc__ = (
        "WARNING: This function is deprecated in favor of minimum_bounding_circle\n"
        + minimum_bounding_circle.__doc__
    )

    def _skyum_lists(points):
        points = points.tolist()
        removed = []
        i = 0
        while True:
            angles = [
                _angle(
                    _prec(p, points),
                    p,
                    _succ(p, points),
                )
                for p in points
            ]
            circles = [
                _circle(
                    _prec(p, points),
                    p,
                    _succ(p, points),
                )
                for p in points
            ]
            radii = [c[0] for c in circles]
            lexord = np.lexsort((radii, angles))  # confusing as hell defaults...
            lexmax = lexord[-1]
            candidate = (
                _prec(points[lexmax], points),
                points[lexmax],
                _succ(points[lexmax], points),
            )
            if angles[lexmax] <= (np.pi / 2.0):
                # print("Constrained by points: {}".format(candidate))
                return _circle(*candidate), points, removed, candidate
            else:
                try:
                    removed.append((points.pop(lexmax), i))
                except IndexError:
                    raise Exception("Construction of Minimum Bounding Circle failed!")
            i += 1

    try:
        from numba import njit, boolean

        HAS_NUMBA = True

        @njit(fastmath=True)
        def _skyum_numba(points):
            i = 0
            complete = False
            while not complete:
                complete, points, candidate, circle = _skyum_iteration(points)
                if complete:
                    return circle, points, None, candidate

        @njit(fastmath=True)
        def _skyum_iteration(points):
            points = points.reshape(-1, 2)
            n = points.shape[0]
            angles = np.empty((n,))
            circles = np.empty((n, 3))
            for i in range(n):
                p = points[(i - 1) % n]
                q = points[i % n]
                r = points[(i + 1) % n]
                angles[i] = _angle(p, q, r)
                circles[i] = _circle(p, q, r)
            radii = circles[:, 0]
            # workaround for no lexsort in numba
            angle_argmax = angles.argmax()
            angle_max = angles[angle_argmax]
            # the maximum radius for the largest angle
            lexmax = (radii * (angles == angle_max)).argmax()

            candidate = (lexmax - 1) % n, lexmax, (lexmax + 1) % n
            if angles[lexmax] <= (np.pi / 2.0):
                return True, points, lexmax, circles[lexmax]
            else:
                mask = np.ones((n,), dtype=boolean)
                mask[lexmax] = False
                new_points = points[mask, :]
                return False, new_points, lexmax, circles[lexmax]

    except ModuleNotFoundError:

        def njit(func, **kwargs):
            return func

    @njit
    def _angle(p, q, r):
        pq = p - q
        rq = r - q
        magnitudes = np.linalg.norm(pq) * np.linalg.norm(rq)
        return np.abs(np.arccos(np.dot(pq, rq) / magnitudes))

    def _prec(p, l):
        """
        retrieve the predecessor of p in list l
        """
        pos = l.index(p)
        if pos - 1 < 0:
            return l[-1]
        else:
            return l[pos - 1]

    def _succ(p, l):
        """
        retrieve the successor of p in list l
        """
        pos = l.index(p)
        if pos + 1 >= len(l):
            return l[0]
        else:
            return l[pos + 1]

    @njit
    def _euclidean_distance(px, py, qx, qy):
        return np.sqrt((px - qx) ** 2 + (py - qy) ** 2)

    @njit
    def _circle(p, q, r, dmetric=_euclidean_distance):
        """
        Returns (radius, (center_x, center_y)) of the circumscribed circle by the
        triangle pqr.

        note, this does not assume that p!=q!=r
        """
        px, py = p
        qx, qy = q
        rx, ry = r
        angle = np.abs(_angle(p, q, r))
        if np.abs(angle - np.pi) < 1e-5:  # angle is pi
            radius = _euclidean_distance(px, py, rx, ry) / 2.0
            center_x = (px + rx) / 2.0
            center_y = (py + ry) / 2.0
        elif np.abs(angle) < 1e-5:  # angle is zero
            radius = _euclidean_distance(px, py, qx, qy) / 2.0
            center_x = (px + qx) / 2.0
            center_y = (py + qy) / 2.0
        else:
            D = 2 * (px * (qy - ry) + qx * (ry - py) + rx * (py - qy))
            center_x = (
                (px**2 + py**2) * (qy - ry)
                + (qx**2 + qy**2) * (ry - py)
                + (rx**2 + ry**2) * (py - qy)
            ) / float(D)
            center_y = (
                (px**2 + py**2) * (rx - qx)
                + (qx**2 + qy**2) * (px - rx)
                + (rx**2 + ry**2) * (qx - px)
            ) / float(D)
            radius = _euclidean_distance(center_x, center_y, px, py)
        return radius, center_x, center_y