# 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

    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

    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 OpenCV2, so
        if that is not available, then this function will fail.

        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).
            Computed directly from cv2.minAreaRect.

        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)
        try:
            from cv2 import minAreaRect, boxPoints
        except (ImportError, ModuleNotFoundError):
            raise ModuleNotFoundError("OpenCV2 is required to use this function. Please install it with `pip install opencv-contrib-python`")
        ((x, y), (w, h), angle) = rot_rect = minAreaRect(points.astype(np.float32))
        out_points = boxPoints(rot_rect)
        if return_angle:
            return (out_points, angle)
        return out_points

    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