# Source code for pypdevsbbl.extra.mathutils

    #  Copyright 2020 Modelling, Simulation and Design Lab (MSDL) at
    #  McGill University and the University of Antwerp (http://msdl.cs.mcgill.ca/)
    #
    #  Licensed under the Apache License, Version 2.0 (the "License");
    #  you may not use this file except in compliance with the License.
    #  You may obtain a copy of the License at
    #
    #     http://www.apache.org/licenses/LICENSE-2.0
    #
    #  Unless required by applicable law or agreed to in writing, software
    #  distributed under the License is distributed on an "AS IS" BASIS,
    #  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    #  See the License for the specific language governing permissions and
    #  limitations under the License.

    """This file contains the main functionality for having more control over boxplots.

    These can be used for custom computation of boxplots, histograms and some other features,
    based upon certain algorithms. This way, it is possible to compute these values in
    situations where you don't have access over the entire dataset, but only over a few values
    thereof.

    There are also some additional helper functions that can be useful in other circumstances
    as well.
    """

    from enum import Enum
    import math

    __all__ = ['BoxplotAlgorithm', 'compute_boxplot', 'FAME', 'histogram', 'histostep', 'PPP', 'linear_interpolation',
               'sign', 'binomial', 'binprob', 'draw_boxplot', 'mean', 'var', 'std', 'harmonic', 'harmonicinv']

    class BoxplotAlgorithm(Enum):
        """Helper enum for creating boxplots of series of :math:`N` items.

        Boxplots are created with 5 values: [`min`, `Q1`, `median`, `Q3`, `max`]

        Even though `Q1` and `Q3` are also defined as being medians, some software
        packages use other algorithms for computing these values. This enum allows
        a quick and clean access of a set of these other algorithms.

        See Also:
            :meth:`linear_interpolation`
        """
        MEDIAN             = 0
        """Default computation where the `Q1` and `Q3` values are the medians of the
        lower list and the upper list of the set split on its global median, respectively.

        The global median is not included in the computation. When there is an even amount
        of items, the two  values that make up the median are both excluded from the
        computation.

        Computes the element at indexes `Q1` and `Q3` and determines a linear interpolation
        on the series if the index is not an integer.

        :math:`Q1 = \\begin{cases}\\dfrac{N}{4}\\textrm{ if }N\\textrm{ is even}\\\\\\\\
        \\dfrac{N + 1}{4}\\textrm{ if }N\\textrm{ is odd}\\end{cases}\\qquad
        Q3 = \\begin{cases}\\dfrac{3\\cdot N + 4}{4}\\textrm{ if }N\\textrm{ is even}\\\\\\\\
        \\dfrac{3\\cdot N + 3}{4}\\textrm{ if }N\\textrm{ is odd}\\end{cases}`
        """

        TUKEY              = 1
        """Computation that's used by `R` in the :meth:`fivenum` function and by
        :mod:`matplotlib` for plotting boxplots.

        The computation is the same as :attr:`MEDIAN`, except that the median is included
        in the computation.

        Computes the element at indexes `Q1` and `Q3` and determines a linear interpolation
        on the series if the index is not an integer.

        :math:`Q1 = \\begin{cases}\\dfrac{N + 2}{4}\\textrm{ if }N\\textrm{ is even}\\\\\\\\
        \\dfrac{N + 3}{4}\\textrm{ if }N\\textrm{ is odd}\\end{cases}\\qquad
        Q3 = \\begin{cases}\\dfrac{3\\cdot N + 2}{4}\\textrm{ if }N\\textrm{ is even}\\\\\\\\
        \\dfrac{3\\cdot N + 1}{4}\\textrm{ if }N\\textrm{ is odd}\\end{cases}`
        """

        MINITAB            = 2
        """Computation that's used by `Minitab` and by `MS Excel` (in any version later
        than 2007, in the :meth:`QUARTILE.EXC` function).

        Computes the element at indexes `Q1` and `Q3` and determines a linear interpolation
        on the series if the index is not an integer.

        :math:`Q1 = \\dfrac{N + 1}{4}\\qquad Q3 = \\dfrac{3\\cdot N + 3}{4}`
        """

        FREUND_PERLES      = 3
        """Computation that's used by `MS Excel` (version 2007 and earlier use :meth:`QUARTILE`,
        which was replaced in later versions by :meth:`QUARTILE.INC`), by `R` in the :meth:`summary`
        function and by :meth:`numpy.percentile` (with a linear interpolation).

        Computes the element at indexes `Q1` and `Q3` and determines a linear interpolation
        on the series if the index is not an integer.

        :math:`Q1 = \\dfrac{N + 3}{4}\\qquad Q3 = \\dfrac{3\\cdot N + 1}{4}`
        """

        MOORE_MCCABE       = 4
        """Similar to :attr:`MEDIAN`, as is it used in the `TI-83` calculators for boxplots.

        Computes the element at indexes `Q1` and `Q3` and determines a linear interpolation
        on the series if the index is not an integer.

        :math:`Q1 = \\begin{cases}\\dfrac{N + 2}{4}\\textrm{ if }N\\textrm{ is even}\\\\\\\\
        \\dfrac{N + 1}{4}\\textrm{ if }N\\textrm{ is odd}\\end{cases}\\qquad
        Q3 = \\begin{cases}\\dfrac{3\\cdot N + 2}{4}\\textrm{ if }N\\textrm{ is even}\\\\\\\\
        \\dfrac{3\\cdot N + 3}{4}\\textrm{ if }N\\textrm{ is odd}\\end{cases}`
        """

        MENDENHALL_SINCICH = 5
        """Same as :attr:`MINITAB`, except this method uses rounding for its computations.

        Computes the element at indexes `Q1` and `Q3` and determines a linear interpolation
        on the series if the index is not an integer.

        :math:`Q1 = \\left\\lfloor\\dfrac{N + 1}{4}\\right\\rfloor\\qquad
        Q3 = \\left\\lceil\\dfrac{3\\cdot N + 3}{4}\\right\\rceil`
        """

        FAME               = 6
        """Computes the `Q1` and `Q3` values via the *Fast Algorithm for Median Estimation*.

        See Also:
            :meth:`FAME`
        """

        P_SQUARE           = 7
        """Computes the `Q1` and `Q3` values via the *P^2 Algorithm for Dynamic Calculation of Quantiles
        and Histograms without Scoring Observations*.

        See Also:
            - :meth:`histogram`
            - :meth:`histostep`
            - :meth:`PPP`
        """

    def compute_boxplot(series: list, algorithm: BoxplotAlgorithm = BoxplotAlgorithm.MEDIAN):
        """Compute the boxplot values for a series.

        This function computes the [`min`, `Q1`, `median`, `Q3`, `max`] set of a series,
        based upon a certain algorithm.

        Note:
             The median in the result is the definition-wise median and not computed with some
             other algorithm.

        See Also:
            :meth:`linear_interpolation`

        Args:
            series (list):  A list of numeric vales from which a boxplot
                            must be computed.

            algorithm (BoxplotAlgorithm):
                            The algorithm to use for computing the data.
                            Even though the `Q1` and `Q3` values are technically defined as
                            medians, some software actually use other computations for these
                            values. Defaults to :attr:`BoxplotAlgorithm.MEDIAN`.

        Returns:
            A dictionary with keys [min, Q1, median, Q3, max] and their corresponding values.
        """
        s = sorted(series)
        N = len(s)
        M = lambda x: (x + 1) / 2      # Compute 1-based median index of series
        res = {
            "min": min(s),
            "max": max(s),
            "median": linear_interpolation(M(N) - 1, s)
        }

        if algorithm is BoxplotAlgorithm.MEDIAN:
            idx1 = (N + (0 if N % 2 == 0 else 1)) / 4
            idx3 = (3 * N + (4 if N % 2 == 0 else 3)) / 4
        elif algorithm is BoxplotAlgorithm.TUKEY:
            idx1 = (N + (2 if N % 2 == 0 else 3)) / 4
            idx3 = (N * 3 + (2 if N % 2 == 0 else 1)) / 4
        elif algorithm is BoxplotAlgorithm.MINITAB:
            idx1 = (N + 1) / 4
            idx3 = (N * 3 + 3) / 4
        elif algorithm is BoxplotAlgorithm.FREUND_PERLES:
            idx1 = (N + 3) / 4
            idx3 = (N * 3 + 1) / 4
        elif algorithm is BoxplotAlgorithm.MOORE_MCCABE:
            idx1 = (N + (2 if N % 2 == 0 else 1)) / 4
            idx3 = (N * 3 + (2 if N % 2 == 0 else 3)) / 4
        elif algorithm is BoxplotAlgorithm.MENDENHALL_SINCICH:
            idx1 = math.floor((N + 1) / 4)
            idx3 = math.ceil((N * 3 + 3) / 4)
        elif algorithm in [BoxplotAlgorithm.FAME, BoxplotAlgorithm.P_SQUARE]:
            # run it afterwards for efficiency, so set dummy data
            # This data will not be used anyways, but it prevents
            #    an error in future lines.
            idx1 = 1
            idx3 = 1
        else:
            raise ValueError("Unknown Algorithm")

        res["Q1"] = linear_interpolation(idx1 - 1, s)
        res["Q3"] = linear_interpolation(idx3 - 1, s)

        if algorithm is BoxplotAlgorithm.FAME:
            Mix = M(N) - 1
            res["Q1"] = FAME(s[:math.floor(Mix)])
            res["Q3"] = FAME(s[math.ceil(Mix)+1:])
        elif algorithm is BoxplotAlgorithm.P_SQUARE:
            x, _ = histogram(series, 4)
            res["Q1"] = x[1]
            res["Q3"] = x[3]

        return res

    def FAME(series: list, b = 0.1):
        """Do the `FAME` algorithm based upon a series.

        Loops throughout the entire series and tries to estimate the median.
        Please note that this is a mere estimate and not at all the actual
        median!

        `FAME` (Fast Algorithm for Median Estimation) is described in
        *An Optimal Median Calculation Algorithm for Estimating Internet Link
        Delays From Active Measurements* by `Dima Feldman` and `Yuval Shavitt`
        at the university of Tel-Aviv, 2007.

        See Also:
            - http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.108.7376&rep=rep1&type=pdf (paper)
            - http://www.eng.tau.ac.il/~shavitt/courses/LargeG/streaming-median.pdf

        Args:
            series (list):  The series to loop through.
            b (numeric):    Minimal initial step size. Defaults to 0.1.

        Returns:
            The estimated median of the series.
        """
        M = series[0]
        step = max(abs(M / 2), b)

        for i in series[1:]:
            if M > i:
                M -= step
            elif M < i:
                M += step
            if abs(i - M) < step:
                step /= 2

        return M

    def histogram(series: list, b: int):
        """Computes a b-cell histogram estimate of the series.

        Loops throughout the series and plots an estimate for
        a b+1 point histogram, using b equiprobable cells.
        When there are less then b observations, the series
        itself becomes the histogram.

        The source for this algorithm is based upon the
        *P^2 Algorithm for Dynamic Calculation of Quantiles
        and Histograms without Scoring Observations* by `Raj
        Jain` and `Imrich Chlamtac` at the Washington University
        in Saint Louis, 1985.

        See Also:
            - https://www.cse.wustl.edu/~jain/papers/ftp/psqr.pdf (paper)
            - :meth:`histostep`
            - :meth:`PPP`

        Args:
            series (list):  The series of which to compute the
                            histogram.
            b (int):        The amount of cells to compute.
                            Note that b+1 will identify the amount
                            of 'bars' in the histogram.

        Returns:
            An (x, y) tuple where x contains the x-locations of all
            values on the x-axis of the histogram and y contains all
            heights of all bars. Both x and y are lists of numbers of
            size b+1.
        """

        if len(series) <= b:
            return sorted(series), [x/len(series) for x in range(len(series))]

        q = sorted(series[:b+1])
        o = series[b+1:]
        n = list(range(b+1))

        for j in range(len(o)):
            xj = o[j]
            N = b + j + 2   # (b+1) observations before loop + currently on (j+1) observation in loop
            histostep(xj, b, q, n, N)

        return q, [x/len(series) for x in n]

    def histostep(xj, b: int, q: list, n: list, N: int):
        """Helper function for the histogram function.

        This function was extracted to make it easier for the algorithm
        to be used without knowledge of the entire series.

        See Also:
            - :meth:`histogram`
            - :meth:`PPP`

        Args:
            xj (numeric):   The new item in the series.
            b (int):        The amount of cells to compute.
            q (list):       The datavalues of the bars in the histogram.
            n (list):       The positions of the bars in the histogram.
            N (int):        The current observation index.
        """
        if xj < q[0]:
            q[0] = xj
            k = 0
        for i in range(b):
            if q[i] <= xj < q[i+1]:
                k = i
                break
        if q[b] < xj:
            q[b] = xj
            k = b-1

        for i in range(k+1, b+1):
            n[i] += 1

        for i in range(1, b):
            # Calculate desired marker position
            n_ = i * (N - 1) / b
            di = n_ - n[i]
            # Prevent markers moving to the same position.
            if (di >= 1 and n[i+1] - n[i] > 1) or \
                (di <= -1 and n[i-1] - n[i] < 1):
                di = sign(di)
                qi = PPP(q, n, i, di)
                if q[i-1] < qi < q[i+1]:
                    q[i] = qi
                else:
                    q[i] += di * (q[i + di] - q[i]) / (n[i + di] - n[i])
                n[i] += di

    def PPP(q: list, n: list, i: int, d=1):
        """Do a Piecewise-Parabolic Prediction.

        Assume we have a parabola passing through
        (n[i-1], q[i-1]), (n[i], q[i]) and (n[i+1], q[i+1])
        which is of the form

        .. math:: y = ax^2 + bx + c

        where (x, y) are coordinates in these points.
        The interpolated value (prediction) at location i can
        be determined via solving the equations set above.

        Based upon the PP-algorithm from
        https://www.cse.wustl.edu/~jain/papers/ftp/psqr.pdf

        Args:
            q (list):   The heights of the points.
            n (list):   The location of the points.
            i (int):    The index where we will look.
            d (int):    Either -1 (left) or 1 (right)
        """
        assert abs(d) == 1

        def int_n(j, f=1):
            return n[j] - n[j-1] + d*f

        def int_q(j):
            return (q[j+1]-q[j]) / (n[j+1] - n[j])

        return q[i] + (d / (n[i+1] - n[i-1])) * (int_n(i) * int_q(i) + int_n(i+1, -1) * int_q(i-1))

    def linear_interpolation(x, data: list):
        """Compute the linear interpolation.

        Using the data, it will compute the value that should appear at index idx.
        This function will do a linear interpolation between the datapoints
        i-1 and i around x. (Whereas i-1=floor(x) and i=ceil(x).)
        We assume the line through these points follows the equation

        .. math::
                    y &= mx + q \\\\
                    y - y1 &= (y2 - y1) / (x2 - x1) * (x - x1)

        Args:
            x (numeric):    The location to do the interpolation at. (the x-value)
            data (list):    The list of y-values at integer values.

        Returns:
            The value that corresponds with the idx in the list.
        """
        if math.floor(x) == x:
            return data[int(x)]

        x0 = math.floor(x)
        y0 = data[x0]
        x1 = math.ceil(x)
        y1 = data[x1]

        return y0 + (x - x0) * (y1 - y0) / (x1 - x0)

    def sign(num):
        """Apply the sign function.

        Args:
            num (numeric):    The value of which the sign needs to be determined.

        Returns:
            :math:`\\begin{cases}\\begin{align}&1&\\textrm{if num > 0}\\\\&-1&\\textrm{if num < 0}\\\\&0&\\textrm{otherwise }\\end{align}\\end{cases}`
        """
        return 1 if num > 0 else -1 if num < 0 else 0

    def binomial(n, k):
        """The binomial computation.

        The amount of possible ways to choose k out of n.
        This is the computation of the binomial coefficient C(n, k).
        Instead of the basic formula of ``n! / (k! * (n - k)!)``
        where 'x!' stands for 'x factorial', this method will try and
        compute a space-efficient binomial (i.e. O(1)).

        It is based on the explanation on:
        https://www.geeksforgeeks.org/space-and-time-efficient-binomial-coefficient/

        Args:
            n (numeric):    The total number (top number).
            k (numeric):    The amount to choose (bottom number).

        Returns:
            The binomial of (n k).
        """
        # We know C(n, k) = C(n, n - k) and therefore we can simplify
        #   Basically, we choose the least amount of operations needed.
        if k > n - k:
            k = n - k
        res = 1.0
        for i in range(k):
            res *= float(n - i)
            res /= float(i + 1)
        return res

    def binprob(t, n, p):
        """Compute the probability w.r.t. a binomial distribution.

        Assume X follows a discrete binomial distribution with size n
        and probability p: X ~ Bin(n, p)

        Args:
            t (numeric):    The value to compute.
            n (numeric):    The size of the distribution.
            p (float):      The probability of the distribution.

        Returns:
            The probability that X is less then t:
                    Pr[X <= t]
        """
        # print("\tX ~ Bin(%f, %f); Pr[X <= %i]" % (n, p, t))
        return sum([binomial(n, i)*(p**i)*((1-p)**(n-i)) for i in range(math.floor(t))])

    def draw_boxplot(minimum, Q1, median, Q3, maximum, outliers=False, labels=None):
        """Get a boxplot-data object in matplotlib, based on its values.

        Using this function, it is possible to quickly draw a boxplot, based on
        its internal values. By default, matplotlib uses either bootstrapping or
        Gaussian-based asymptotic approximation (Tukey).

        This function therefore allows drawing boxplots based on other algorithms,
        but it does NOT draw the boxplots themselves. It yields a mere object, which
        can be used in the Axes.bpx method of matplotlib.

        Args:
            minimum (numeric):  The smallest number to be drawn on the boxplot.
            Q1 (numeric):       The lower quartile edge of the box.
            median (numeric):   The median position indicator in the box.
            Q3 (numeric):       The upper quartile edge of the box.
            maximum (numeric):  The largest number to be drawn on the boxplot.
            outliers (bool):    When true, this will make sure the whiskers are at most
                                1.5 * IQR away from the box. If the `minimum` and the
                                `maximum` are outside of this range, they will be
                                indicated as outliers.
                                Do note that, as was the case with the compute_boxplot
                                function, this method is a mere approximation of the
                                boxplot, which becomes clear when comparing it to the
                                actual boxplot. Theoretically, it's a good approximation
                                nonetheless.
                                Defaults to ``False``.
            labels (iterable):  Labels for each dataset.
                                See also cbook.boxplot_stats.labels

        Returns:
            A dictionary containing the results for each column of data.
            This is basically a transformed dictionary, based on:

            >>> cbook.boxplot_stats([minimum, Q1, median, Q3, maximum], labels=labels)[0]

        Example:
            >>> import matplotlib.pyplot as plt
            >>>
            >>> fig, ax = plt.subplots()
            >>> boxplot = draw_boxplot(minimum, Q1, median, Q3, maximum)
            >>> ax.bpx(boxplot)       # draws the actual boxplot
            >>> plt.show()
        """
        from matplotlib import cbook

        s = cbook.boxplot_stats([minimum, Q1, median, Q3, maximum], labels=labels)[0]
        iqr = Q3 - Q1
        s["q1"] = Q1
        s["median"] = median
        s["q3"] = Q3
        s["iqr"] = iqr
        s["fliers"] = []
        if outliers:
            L = Q1 - (1.5 * iqr)
            H = Q3 + (1.5 * iqr)
            if minimum < L:
                s["fliers"].append(minimum)
            if maximum > H:
                s["fliers"].append(maximum)
            s["whislo"] = max(minimum, L)
            s["whishi"] = min(maximum, H)
        else:
            s["whislo"] = minimum
            s["whishi"] = maximum

        return [s]

    def mean(seq):
        """Computes the mean of a sequence.

        Args:
            seq (iterable): The sequence to compute the mean for.

        Returns:
            The mean of the sequence.
        """
        if len(seq) == 0:
            return 0
        return sum(seq) / len(seq)

    def var(seq):
        """Computes the variance of a sequence.

        Args:
            seq (iterable): The sequence to compute the variance for.

        Returns:
            The variance of the sequence.
        """
        mu = mean(seq)
        S = sum([(X - mu)**2 for X in seq])
        return S / len(seq)

    def std(seq):
        """Computes the standard derivation of a sequence.

        Args:
            seq (iterable): The sequence to compute the standard
                            derivation for.

        Returns:
            The standard derivation of the sequence.
        """
        return var(seq) ** 0.5

    def harmonic(n: int, m=1):
        """Computes the n-th generalized harmonic number.

        Args:
            n (int):        The number to obtain.
            m (numeric):    The power of each denominator.
                            Defaults to 1, which also allows non-generalized numbers.

        Returns:
            The n-th generalized harmonic number.
        """
        return sum([1 / (q ** m) for q in range(1, n+1)])

    def harmonicinv(y, m=1):
        """Compute the inverse general harmonic number.

        I.e. :math:`GH(n, m) = y \\iff IH(y, m) = n`

        Args:
            y (numeric):    The general harmonic number.
            m (numeric):    The power to compute.
                            Defaults to 1.

        Returns:
            The inverse general harmonic number.
        """
        s = 0
        while y > 0:
            s += 1
            y -= 1 / (s ** m)
        return s

    if __name__ == '__main__':
        # l = [0, 0, 1, 2, 63, 61, 27, 13]
        # l = [0, 1, 2, 4, 7, 8, 9, 10, 12]
        l = list(range(1, 102))
        algos = [
            BoxplotAlgorithm.MEDIAN,
            BoxplotAlgorithm.TUKEY,
            BoxplotAlgorithm.MINITAB,
            BoxplotAlgorithm.FREUND_PERLES,
            BoxplotAlgorithm.MOORE_MCCABE,
            BoxplotAlgorithm.MENDENHALL_SINCICH,
            BoxplotAlgorithm.FAME,
            BoxplotAlgorithm.P_SQUARE,
        ]
        for algo in algos:
            print(algo.name, compute_boxplot(l, algo))