def ravel_indices(indices, shape):
    """
    Convert nD to 1D indices for an array of given shape.
        flat_indices = ravel_indices(indices, size)
    
    :Input:
        indices: array of indices. Should be integer and have shape=([S],D), 
                 for S the "subshape" of indices array, pointing to a D dimensional array.
        shape:   shape of the nd-array these indices are pointing to (a tuple/list/ of length D)
        
    :Output: 
        flat_indices: an array of shape S
    
    :Note: 
       This is the opposite of unravel_indices: for any tuple 'shape'
          ind is equal to    ravel_indices(unravel_indices(ind,shape),shape)
                   and to  unravel_indices(  ravel_indices(ind,shape),shape)
    """
    dim_prod = _np.cumprod([1] + list(shape)[:0:-1])[_np.newaxis,::-1]
    ind = _np.asanyarray(indices)
    S   = ind.shape[:-1]
    K   = _np.asanyarray(shape).size
    ind.shape = S + (K,)
    
    return _np.sum(ind*dim_prod,-1)

    def fit(self, X, y, **params):
        """
        Fit Ridge regression model

        Parameters
        ----------
        X : numpy array of shape [n_samples,n_features]
            Training data
        y : numpy array of shape [n_samples]
            Target values

        Returns
        -------
        self : returns an instance of self.
        """
        self._set_params(**params)

        X = np.asanyarray(X, dtype=np.float)
        y = np.asanyarray(y, dtype=np.float)

        n_samples, n_features = X.shape

        X, y, Xmean, ymean = self._center_data(X, y)

        if n_samples > n_features:
            # w = inv(X^t X + alpha*Id) * X.T y
            self.coef_ = linalg.solve(np.dot(X.T, X) + self.alpha * np.eye(n_features), np.dot(X.T, y))
        else:
            # w = X.T * inv(X X^t + alpha*Id) y
            self.coef_ = np.dot(X.T, linalg.solve(np.dot(X, X.T) + self.alpha * np.eye(n_samples), y))

        self._set_intercept(Xmean, ymean)
        return self

 def eval(self, x, y):
     """Evaluate model at a given position ``(x, y)`` position.
     """
     x = np.asanyarray(x, dtype=float)
     y = np.asanyarray(y, dtype=float)
     parvals = self.parvals(x)
     return self._eval_y(y, parvals)

    def __call__(self, projectables, **info):
        if len(projectables) != 2:
            raise ValueError("Expected 2 datasets, got %d" %
                             (len(projectables), ))

        # TODO: support datasets with palette to delegate this to the image
        # writer.

        data, palette = projectables
        palette = np.asanyarray(palette).squeeze() / 255.0
        colormap = self.build_colormap(palette, data.dtype, data.attrs)

        channels, colors = colormap.palettize(np.asanyarray(data.squeeze()))
        channels = palette[channels]
        fill_value = data.attrs.get('_FillValue', np.nan)
        if np.isnan(fill_value):
            mask = data.notnull()
        else:
            mask = data != data.attrs['_FillValue']
        r = xr.DataArray(channels[:, :, 0].reshape(data.shape),
                         dims=data.dims, coords=data.coords,
                         attrs=data.attrs).where(mask)
        g = xr.DataArray(channels[:, :, 1].reshape(data.shape),
                         dims=data.dims, coords=data.coords,
                         attrs=data.attrs).where(mask)
        b = xr.DataArray(channels[:, :, 2].reshape(data.shape),
                         dims=data.dims, coords=data.coords,
                         attrs=data.attrs).where(mask)

        res = super(PaletteCompositor, self).__call__((r, g, b), **data.attrs)
        res.attrs['_FillValue'] = np.nan
        return res

    def fit(self, X, y):
        """
        Fit linear model.

        Parameters
        ----------
        X : numpy array of shape [n_samples,n_features]
            Training data
        y : numpy array of shape [n_samples]
            Target values
        fit_intercept : boolean, optional
            wether to calculate the intercept for this model. If set
            to false, no intercept will be used in calculations
            (e.g. data is expected to be already centered).
        normalize : boolean, optional
            If True, the regressors X are normalized

        Returns
        -------
        self : returns an instance of self.
        """
        X = np.asanyarray(X)
        y = np.asanyarray(y)

        X, y, X_mean, y_mean, X_std = self._center_data(X, y,
                self.fit_intercept, self.normalize, self.copy_X)

        self.coef_, self.residues_, self.rank_, self.singular_ = \
                np.linalg.lstsq(X, y)

        self._set_intercept(X_mean, y_mean, X_std)
        return self

    def __init__(self, points, inside=None):
        r"""
        Parameters
        ----------
        points : An Nx3 array of (*x*, *y*, *z*) triples in vector space
            These points define the boundary of the polygon.  It must
            be "closed", i.e., the last point is the same as the first.

            It may contain zero points, in which it defines the null
            polygon.  It may not contain one, two or three points.
            Four points are needed to define a triangle, since the
            polygon must be closed.

        inside : An (*x*, *y*, *z*) triple, optional
            This point must be inside the polygon.  If not provided, the
            mean of the points will be used.
        """
        if len(points) == 0:
            # handle special case of initializing with an empty list of
            # vertices (ticket #1079).
            self._inside = np.zeros(3)
            self._points = np.asanyarray(points)
            return
        elif len(points) < 3:
            raise ValueError("Polygon made of too few points")
        else:
            assert np.array_equal(points[0], points[-1]), 'Polygon is not closed'

        self._points = points = np.asanyarray(points)

        if inside is None:
            self._inside = self._find_new_inside(points)
        else:
            self._inside = np.asanyarray(inside)

def chi2(N_S, B, S, sigma2):
    r"""Chi-square statistic with user-specified variance.

     .. math::

         \chi^2 = \frac{(N_S - B - S) ^ 2}{\sigma ^ 2}

    Parameters
    ----------
    N_S : array_like
        Number of observed counts
    B : array_like
        Model background
    S : array_like
        Model signal
    sigma2 : array_like
        Variance

    Returns
    -------
    stat : ndarray
        Statistic per bin

    References
    ----------
    * Sherpa stats page (http://cxc.cfa.harvard.edu/sherpa/statistics/#chisq)
    """
    N_S = np.asanyarray(N_S, dtype=np.float64)
    B = np.asanyarray(B, dtype=np.float64)
    S = np.asanyarray(S, dtype=np.float64)
    sigma2 = np.asanyarray(sigma2, dtype=np.float64)

    stat = (N_S - B - S) ** 2 / sigma2

    return stat

    def fit(self, X, y, **params):
        """
        Fit linear model.

        Parameters
        ----------
        X : numpy array of shape [n_samples,n_features]
            Training data
        y : numpy array of shape [n_samples]
            Target values
        fit_intercept : boolean, optional
            wether to calculate the intercept for this model. If set
            to false, no intercept will be used in calculations
            (e.g. data is expected to be already centered).

        Returns
        -------
        self : returns an instance of self.
        """
        self._set_params(**params)
        X = np.asanyarray(X)
        y = np.asanyarray(y)

        X, y, Xmean, ymean = LinearModel._center_data(X, y, self.fit_intercept)

        self.coef_, self.residues_, self.rank_, self.singular_ = \
                np.linalg.lstsq(X, y)

        self._set_intercept(Xmean, ymean)
        return self

    def predict(self, X, copy=True):
        """Apply the dimension reduction learned on the train data.
            Parameters
            ----------
            X: array-like of predictors, shape (n_samples, p)
                Training vectors, where n_samples in the number of samples and
                p is the number of predictors.

            copy: X has to be normalize, do it on a copy or in place
                with side effect!

            Notes
            -----
            This call require the estimation of a p x q matrix, which may
            be an issue in high dimensional space.
        """
        # Normalize
        if copy:
            Xc = (np.asanyarray(X) - self.x_mean_)
        else:
            X = np.asanyarray(X)
            Xc -= self.x_mean_
            Xc /= self.x_std_
        Ypred = np.dot(Xc, self.coefs)
        return Ypred + self.y_mean_

def background_error(n_off, alpha):
    r"""Estimate standard error on background
    in the on region from an off-region observation.

    .. math::

          \Delta\mu_{bkg} = \alpha \times \sqrt{n_{off}}

    Parameters
    ----------
    n_off : array_like
        Observed number of counts in the off region
    alpha : array_like
        On / off region exposure ratio for background events

    Returns
    -------
    background : ndarray
        Background estimate for the on region

    Examples
    --------
    >>> background_error(n_off=4, alpha=0.1)
    0.2
    >>> background_error(n_off=9, alpha=0.2)
    0.6
    """
    n_off = np.asanyarray(n_off, dtype=np.float64)
    alpha = np.asanyarray(alpha, dtype=np.float64)

    return alpha * sqrt(n_off)

def excess(n_on, n_off, alpha):
    r"""Estimate excess in the on region for an on-off observation.

    .. math::

          \mu_{excess} = n_{on} - \alpha \times n_{off}

    Parameters
    ----------
    n_on : array_like
        Observed number of counts in the on region
    n_off : array_like
        Observed number of counts in the off region
    alpha : array_like
        On / off region exposure ratio for background events

    Returns
    -------
    excess : ndarray
        Excess estimate for the on region

    Examples
    --------
    >>> excess(n_on=10, n_off=20, alpha=0.1)
    8.0
    >>> excess(n_on=4, n_off=9, alpha=0.5)
    -0.5
    """
    n_on = np.asanyarray(n_on, dtype=np.float64)
    n_off = np.asanyarray(n_off, dtype=np.float64)
    alpha = np.asanyarray(alpha, dtype=np.float64)

    return n_on - alpha * n_off

def background(n_off, alpha):
    r"""Estimate background in the on-region from an off-region observation.

    .. math::

        \mu_{background} = \alpha \times n_{off}

    Parameters
    ----------
    n_off : array_like
        Observed number of counts in the off region
    alpha : array_like
        On / off region exposure ratio for background events

    Returns
    -------
    background : ndarray
        Background estimate for the on region

    Examples
    --------
    >>> background(n_off=4, alpha=0.1)
    0.4
    >>> background(n_off=9, alpha=0.2)
    1.8
    """
    n_off = np.asanyarray(n_off, dtype=np.float64)
    alpha = np.asanyarray(alpha, dtype=np.float64)

    return alpha * n_off

    def fit(self, X, y, **params):
        """Fit the model using X, y as training data

        Parameters
        ----------
        X : array-like, shape = [n_samples, n_features]
            Training data.

        y : array-like, shape = [n_samples]
            Target values, array of integer values.

        params : list of keyword, optional
            Overwrite keywords from __init__
        """
        X = np.asanyarray(X)
        if y is None:
            raise ValueError("y must not be None")
        self._y = np.asanyarray(y)
        self._set_params(**params)

        if self.algorithm == 'ball_tree' or \
           (self.algorithm == 'auto' and X.shape[1] < 20):
            self.ball_tree = BallTree(X, self.window_size)
        else:
            self.ball_tree = None
            self._fit_X = X
        return self

def unravel_indices(indices,shape):
    """
    Convert indices in a flatten array to nD indices of the array with given shape.
        nd_indices = unravel_indices(indices, shape)
    
    :Input:
        indices: array/list/tuple of flat indices. Should be integer, of any shape S
        shape:   nD shape of the array these indices are pointing to
        
    :Output: 
        nd_indices: a nd-array of shape [S]xK, where 
                    [S] is the shape of indices input argument
                    and K the size (number of element) of shape     
    
    :Note:
        The algorithm has been inspired from numpy.unravel_index 
        and can be seen as a generalization that manage set of indices
        However, it does not return tuples and no assertion is done on 
        the input indices before convertion:
        The output indices might be negative or bigger than the array size
        
        This is the opposite of ravel_indices:  for any tuple 'shape'
          ind is equal to    ravel_indices(unravel_indices(ind,shape),shape)
                   and to  unravel_indices(  ravel_indices(ind,shape),shape)
    """

    dim_prod = _np.cumprod([1] + list(shape)[:0:-1])[::-1]
    ind = _np.asanyarray(indices)
    S   = ind.shape
    K   = _np.asanyarray(shape).size
    
    ndInd = ind.ravel()[:,_np.newaxis]/dim_prod % shape
    ndInd.shape = S + (K,)
    return ndInd

def interp(x, xp, fp, left=None, right=None):
    """
    Like numpy.interp except for handling of running sequences of
    same values in `xp`.
    """
    x = numpy.asanyarray(x)
    xp = numpy.asanyarray(xp)
    fp = numpy.asanyarray(fp)

    if xp.shape != fp.shape:
        raise ValueError("xp and fp must have the same shape")

    ind = numpy.searchsorted(xp, x, side="right")
    fx = numpy.zeros(len(x))

    under = ind == 0
    over = ind == len(xp)
    between = ~under & ~over

    fx[under] = left if left is not None else fp[0]
    fx[over] = right if right is not None else fp[-1]

    if right is not None:
        # Fix points exactly on the right boundary.
        fx[x == xp[-1]] = fp[-1]

    ind = ind[between]

    df = (fp[ind] - fp[ind - 1]) / (xp[ind] - xp[ind - 1])

    fx[between] = df * (x[between] - xp[ind]) + fp[ind]

    return fx

    def fit(self, X, y):
        X = np.asanyarray(X, dtype='d')
        y = np.asanyarray(y, dtype='d')
        
        n = X.shape[0]
        num_dists = self.num_dists
        
        if self.num_dists > n:
            raise ParameterException('Number of distances is greater than ' + \
                    'num rows in X')

        if self.num_dists <= 0:
            self.R = None
        else:
            rand_idx = np.random.choice(X.shape[0], int(num_dists), replace=False)
            self.R = X[rand_idx]
            
            D = np.exp(-1.0 * ((cdist(X, self.R) ** 2) / (2 * (self.sigma ** 2))))
            X = np.hstack((X, D))

        #Un-comment for mrse code
        #X, y = mrse_transform(X, y)

        self.model = self.base_learner.fit(X, y)
        return self

    def get_loss(self, state, action, reward, state_prime, episode_end):
        s = Variable(cuda.to_gpu(state))
        s_dash = Variable(cuda.to_gpu(state_prime))

        q = self.model.q_function(s)  # Get Q-value

        # Generate Target Signals
        tmp = self.model_target.q_function(s_dash)  # Q(s',*)
        tmp = list(map(np.max, tmp.data))  # max_a Q(s',a)
        max_q_prime = np.asanyarray(tmp, dtype=np.float32)
        target = np.asanyarray(copy.deepcopy(q.data.get()), dtype=np.float32)

        for i in range(self.replay_size):
            if episode_end[i][0] is True:
                tmp_ = np.sign(reward[i])
            else:
                #  The sign of reward is used as the reward of DQN!
                tmp_ = np.sign(reward[i]) + self.gamma * max_q_prime[i]

            target[i, action[i]] = tmp_

        # TD-error clipping
        td = Variable(cuda.to_gpu(target)) - q  # TD error
        td_tmp = td.data + 1000.0 * (abs(td.data) <= 1)  # Avoid zero division
        td_clip = td * (abs(td.data) <= 1) + td/abs(td_tmp) * (abs(td.data) > 1)

        zero_val = Variable(cuda.to_gpu(np.zeros((self.replay_size, self.n_act), dtype=np.float32)))
        loss = F.mean_squared_error(td_clip, zero_val)
        return loss, q

def rv2coe(k, r, v):
    """Converts r, v to classical orbital elements.

    This is a wrapper around rv2coe from ast2body.for.

    Parameters
    ----------
    k : float
        Standard gravitational parameter (km^3 / s^2).
    r : array
        Position vector (km).
    v : array
        Velocity vector (km / s).

    Examples
    --------
    Vallado 2001, example 2-5
    >>> r = [6524.834, 6862.875, 6448.296]
    >>> v = [4.901327, 5.533756, -1.976341]
    >>> k = 3.986e5
    >>> rv2coe(k, r, v)
    (36127.55012131963, 0.83285427644495158, 1.5336055626394494,
    3.9775750028016947, 0.93174413995595795, 1.6115511711293014)

    """
    # TODO: Extend for additional arguments arglat, truelon, lonper
    r = np.asanyarray(r).astype(np.float)
    v = np.asanyarray(v).astype(np.float)
    _, a, ecc, inc, omega, argp, nu, _, _, _, _ = _ast2body.rv2coe(r, v, k)
    coe = np.vstack((a, ecc, inc, omega, argp, nu))
    if coe.shape[-1] == 1:
        coe = coe[:, 0]
    return coe

    def set_value(self, value, force=False):
        # Record new value and increment counter

        # Value can't be updated if observed=True
        if self.observed and not force:
            raise AttributeError('Stochastic '+self.__name__+'\'s value cannot be updated if observed flag is set')

        if self.verbose > 0:
            print_('\t' + self.__name__ + ': value set to ', value)

        # Save current value as last_value
        # Don't copy because caching depends on the object's reference.
        self.last_value = self._value

        if self.mask is None:

            if self.dtype.kind != 'O':
                self._value = asanyarray(value, dtype=self.dtype)
                self._value.flags['W']=False
            else:
                self._value = value

        else:

            new_value = self.value.copy()

            new_value[self.mask] = asanyarray(value, dtype=self.dtype)[self.mask]
            self._value = new_value

        self.counter.click()

def kepler(k, r0, v0, tof):
    """Propagates orbit.

    This is a wrapper around kepler from ast2body.for.

    Parameters
    ----------
    k : float
        Gravitational constant of main attractor (km^3 / s^2).
    r0 : array
        Initial position (km).
    v0 : array
        Initial velocity (km).
    tof : float
        Time of flight (s).

    Raises
    ------
    RuntimeError
        If the status of the subroutine is not 'ok'.

    """
    r0 = np.asanyarray(r0).astype(np.float)
    v0 = np.asanyarray(v0).astype(np.float)
    tof = float(tof)
    assert r0.shape == (3,)
    assert v0.shape == (3,)
    r, v, error = _ast2body.kepler(r0, v0, tof, k)
    error = error.strip().decode('ascii')
    if error != 'ok':
        raise RuntimeError("There was an error: {}".format(error))
    return r, v