# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
''' General matrix and other utilities for statistics '''
from __future__ import print_function
from __future__ import absolute_import

import numpy as np
import scipy.interpolate

__docformat__ = 'restructuredtext'


def mad(a, c=0.6745, axis=0):
    """
    Median Absolute Deviation:

    median(abs(a - median(a))) / c
    """
    _shape = a.shape
    a.shape = np.product(a.shape, axis=0)
    m = np.median(np.fabs(a - np.median(a))) / c
    a.shape = _shape
    return m


class StepFunction(object):
    """
    A basic step function: values at the ends are handled in the simplest
    way possible: everything to the left of x[0] is set to ival; everything
    to the right of x[-1] is set to y[-1].

    Examples
    --------
    >>> x = np.arange(20)
    >>> y = np.arange(20)
    >>> f = StepFunction(x, y)
    >>>
    >>> print(f(3.2))
    3.0
    >>> print(f([[3.2,4.5],[24,-3.1]]))
    [[  3.   4.]
     [ 19.   0.]]
    """

    def __init__(self, x, y, ival=0., sorted=False):

        _x = np.asarray(x)
        _y = np.asarray(y)

        if _x.shape != _y.shape:
            raise ValueError(
                'in StepFunction: x and y do not have the same shape')
        if len(_x.shape) != 1:
            raise ValueError('in StepFunction: x and y must be 1-dimensional')

        self.x = np.hstack([[- np.inf], _x])
        self.y = np.hstack([[ival], _y])

        if not sorted:
            asort = np.argsort(self.x)
            self.x = np.take(self.x, asort, 0)
            self.y = np.take(self.y, asort, 0)
        self.n = self.x.shape[0]

    def __call__(self, time):
        tind = np.searchsorted(self.x, time) - 1
        return self.y[tind]


def ECDF(values):
    """
    Return the ECDF of an array as a step function.
    """
    x = np.array(values, copy=True)
    x.sort()
    x.shape = np.product(x.shape, axis=0)
    n = x.shape[0]
    y = (np.arange(n) + 1.) / n
    return StepFunction(x, y)


def monotone_fn_inverter(fn, x, vectorized=True, **keywords):
    """
    Given a monotone function x (no checking is done to verify monotonicity)
    and a set of x values, return an linearly interpolated approximation
    to its inverse from its values on x.
    """

    if vectorized:
        y = fn(x, **keywords)
    else:
        y = []
        for _x in x:
            y.append(fn(_x, **keywords))
        y = np.array(y)

    a = np.argsort(y)

    return scipy.interpolate.interp1d(y[a], x[a])
