r"""
.. _sec-costl1:

Least absolute deviation
====================================================================================================

Description
----------------------------------------------------------------------------------------------------

This cost function detects changes in the median of a signal.
Overall, it is a robust estimator of a shift in the central point (mean, median, mode) of a distribution :cite:`c1-Bai1995`.
Formally, for a signal :math:`\{y_t\}_t` on an interval :math:`I`,

    .. math:: c(y_{I}) = \sum_{t\in I} \|y_t - \bar{y}\|_1

where :math:`\bar{y}` is the componentwise median of :math:`\{y_t\}_{t\in I}`.

Usage
----------------------------------------------------------------------------------------------------

Start with the usual imports and create a signal.

.. code-block:: python

    import numpy as np
    import matplotlib.pylab as plt
    import ruptures as rpt
    # creation of data
    n, dim = 500, 3  # number of samples, dimension
    n_bkps, sigma = 3, 5  # number of change points, noise standart deviation
    signal, bkps = rpt.pw_constant(n, dim, n_bkps, noise_std=sigma)

Then create a :class:`CostL1` instance and print the cost of the sub-signal :code:`signal[50:150]`.

.. code-block:: python

    c = rpt.costs.CostL1().fit(signal)
    print(c.error(50, 150))


You can also compute the sum of costs for a given list of change points.

.. code-block:: python

    print(c.sum_of_costs(bkps))
    print(c.sum_of_costs([10, 100, 200, 250, n]))


In order to use this cost class in a change point detection algorithm (inheriting from :class:`BaseEstimator`), either pass a :class:`CostL1` instance (through the argument ``'custom_cost'``) or set :code:`model="l1"`.

.. code-block:: python

    c = rpt.costs.CostL1(); algo = rpt.Dynp(custom_cost=c)
    # is equivalent to
    algo = rpt.Dynp(model="l1")


Code explanation
----------------------------------------------------------------------------------------------------

.. autoclass:: ruptures.costs.CostL1
    :members:
    :special-members: __init__


.. rubric:: References

.. bibliography:: ../biblio.bib
    :style: alpha
    :cited:
    :labelprefix: C1
    :keyprefix: c1-

"""
import numpy as np
# from scipy.stats import rankdata
from scipy.stats.mstats import rankdata

from ruptures.base import BaseCost
from ruptures.costs import NotEnoughPoints


class CostRank(BaseCost):
    r"""
    Least absolute deviation.
    """

    model = "rank"

    def __init__(self):
        self.sigma = None
        self.ranks = None

    def fit(self, signal):
        """Set parameters of the instance.

        Args:
            signal (array): signal. Shape (n_samples,) or (n_samples, n_features)

        Returns:
            self
        """
        def func(row):
            row = np.reshape((row + 0.5), (-1, 1))
            return row @ row.T

        ranks = np.argsort(signal, axis=0)
        centered_ranks = (ranks - len(signal) / 2).astype(int)

        sigma = np.apply_along_axis(
            func,
            0,
            centered_ranks
        )
        sigma = np.sum(sigma, axis=1)

        self.sigma = sigma
        self.ranks = centered_ranks

        return self

    def error(self, start, end):
        """Return the approximation cost on the segment [start:end].

        Args:
            start (int): start of the segment
            end (int): end of the segment

        Returns:
            float: segment cost

        Raises:
            NotEnoughPoints: when the segment is too short (less than ``'min_size'`` samples).
        """
        mean = np.mean(self.ranks[start:end])
        return -(end - start) * mean.T * self.sigma[start:end] * mean

@@ -2305,16 +2305,45 @@
 rankdata
+%0Afrom numpy.linalg import inv
 %0A%0Afrom r
@@ -2585,16 +2585,42 @@
 s = None
+%0A        self.min_size = 2
 %0A%0A    de
@@ -2819,32 +2819,105 @@
 elf%0A        %22%22%22%0A
+        if signal.ndim == 1:%0A            signal = signal.reshape(-1, 1)%0A%0A
         def func
@@ -3180,17 +3180,17 @@
         
-0
+1
 ,%0A      
@@ -3255,17 +3255,17 @@
 a, axis=
-1
+0
 )%0A%0A     
@@ -3715,32 +3715,107 @@
 s).%0A        %22%22%22%0A
+        if end - start %3C self.min_size:%0A            raise NotEnoughPoints%0A%0A
         mean = n
@@ -3813,16 +3813,27 @@
  mean = 
+np.reshape(
 np.mean(
@@ -3853,16 +3853,34 @@
 art:end%5D
+, axis=0), (-1, 1)
 )%0A      
@@ -3912,18 +3912,22 @@
  mean.T 
-* 
+@ inv(
 self.sig
@@ -3932,23 +3932,31 @@
 igma
-%5Bstart:end%5D * mean
+) @ mean / len(self.ranks)
 %0A

Fix matrix shapes