def test_rdmd_complex64(self):
     m, n = 9, 7
     a = np.array(np.fliplr(np.vander(np.random.rand(m)+1, n)) + 1j*np.fliplr(np.vander(np.random.rand(m)+1, n)), 
                  np.complex64, order='F')
     a_gpu = gpuarray.to_gpu(a)
     f_gpu, b_gpu, v_gpu = rlinalg.rdmd(a_gpu, k=(n-1), p=0, q=1, modes='standard')
     assert np.allclose(a[:,:(n-1)], np.dot(f_gpu.get(), np.dot(np.diag(b_gpu.get()), v_gpu.get()) ), atol_float32)

def plot_results(x, y, yerr, samples, truth=True, color="r", data_fig=None,
                 show=True):
    if data_fig is None:
        # Plot the data.
        data_fig = plot_data(x, y, yerr, truth=truth)
        data_fig, data_ax = _get_fig_ax(data_fig)
    else:
        data_ax = data_fig.gca()

    # Generate the constraints in data space.
    x0 = np.linspace(-5, 5, 500)
    samples = np.atleast_1d(samples)
    if len(samples.shape) == 2:
        lines = np.dot(np.vander(x0, 2), samples[:, :2].T)
        q = np.percentile(lines, [16, 84, 2.5, 97.5], axis=1)
        data_ax.fill_between(x0, q[2], q[3], color=color, alpha=0.1)
        data_ax.fill_between(x0, q[0], q[1], color=color, alpha=0.3)
    else:
        data_ax.plot(x0, np.dot(np.vander(x0, 2), samples[:2]), color=color)

    if show:
        # Plot the triangle plot.
        true = load_data("line_true_params.txt")
        true[2:] = np.log(true[2:])
        triangle_fig = triangle.corner(samples, bins=24,
                                       labels=["m", "b", "ln(alpha)",
                                               "ln(ell)"],
                                       truths=true)
    else:
        triangle_fig = None

    _format_axes(data_ax)
    return data_fig, triangle_fig

def polyfit2dPure(x, y, z, order=2, w=None):
    '''
    References:
    http://pingswept.org/2009/06/15/least-squares-fit-of-a-surface-to-a-3d-cloud-of-points-in-python-(with-ridiculous-application)/
    '''
#    x = np.asarray(x) + 0.0
#    y = np.asarray(y) + 0.0
#    z = np.asarray(z) + 0.0

    deg = order + 1
    Gx = np.vander(x, deg)
    Gy = np.vander(y, deg)
    G = np.hstack((Gx, Gy))

    del x, y, Gx, Gy
    
    # Apply weighting
    if w is not None:
#        w = np.asarray(w) + 0.0
        G *= w[:, np.newaxis]
        z *= w
        
    del w

    m, _, _, _ = np.linalg.lstsq(G, z)
    return m

def noisy_quad_fit(order, Lambda, n_train=20, n_test=81):
    """
      Creates n_train training data points with noise, fits to poly of order,
      then tests on n_test points (noise free).  Uses offset quadratic
    """
    low_x = -2;
    high_x = 2;
    plt.close('all');
    
    train_x = np.linspace(low_x, high_x, n_train);
    X = np.vander(train_x, N = order+1);
    y = (1+train_x**2) + 0.6*(np.random.rand(n_train) - 0.5);
    #y = (np.sin(3*train_x) - (train_x * np.cos(2*train_x))) + 0.6*(np.random.rand(n_train) - 0.5);
    #y = (np.sin(3*train_x) - (train_x * np.cos(2*train_x)));
    #y = (1+train_x**2);
    theta = regress_theta(X,y,Lambda);
    predict_y = np.dot(X,theta);
    print 'Training Error = ', np.max(np.abs(y - predict_y));
    #trainingerror = np.max(np.abs(y - predict_y));
    
    test_x = np.linspace(low_x, high_x, n_test);
    Xt = np.vander(test_x, N = order+1);
    yt = 1+test_x**2;
    #yt = np.sin(3*test_x) - (test_x * np.cos(2*test_x));
    predict_yt = np.dot(Xt,theta);
    print 'Testing Error = ', np.max(np.abs(yt - predict_yt));
    #testingerror = np.max(np.abs(yt - predict_yt));
    
    plt.plot(train_x, y, 'ro');
    plt.plot(train_x, predict_y, 'rx');
    plt.plot(test_x, predict_yt, 'bx');
    plt.show();

 def test_dmd_complex128(self):
     m, n = 9, 7
     a = np.array(np.fliplr(np.vander(np.random.rand(m)+1, n)) + 1j*np.fliplr(np.vander(np.random.rand(m), n)), 
                  np.complex128, order='F')
     a_gpu = gpuarray.to_gpu(a)
     f_gpu, b_gpu, v_gpu, omega = linalg.dmd(a_gpu, modes='standard', return_amplitudes=True, return_vandermonde=True)
     assert np.allclose(a[:,:(n-1)], np.dot(f_gpu.get(), np.dot(np.diag(b_gpu.get()), v_gpu.get()) ), atol_float64)

    def _regressor(self, x):
        """Form normalised regressor / design matrix from set of input vectors.

        Parameters
        ----------
        x : array of float, shape (2, N)
            Input to function as a 2-D numpy array

        Returns
        -------
        X : array of float, shape (P, N)
            Regressor / design matrix to be used in least-squares fit

        Notes
        -----
        This normalises the 2-D input vectors by centering and scaling them.
        It then forms a regressor matrix with a column per input vector. Each
        column is given by the outer product of the monomials of the first
        dimension with the monomials of the second dimension of the input vector,
        in decreasing polynomial order. For example, if *degrees* is (1, 2) and
        the normalised elements of each input vector in *x* are *x_0* and *x_1*,
        respectively, the column takes the form::

            outer([x_0, 1], [x1 ^ 2, x1, 1])
            = [x_0 * x_1 ^ 2, x_0 * x_1, x_0 * 1, 1 * x_1 ^ 2, 1 * x_1, 1 * 1]
            = [x_0 * x_1 ^ 2, x_0 * x_1, x_0, x_1 ^ 2, x_1, 1]

        This is closely related to the Vandermonde matrix of *x*.

        """
        x_norm = (x - self._mean[:, np.newaxis]) / self._scale[:, np.newaxis]
        v1 = np.vander(x_norm[0], self.degrees[0] + 1)
        v2 = np.vander(x_norm[1], self.degrees[1] + 1).T
        return np.vstack([v1[:, n][np.newaxis, :] * v2 for n in xrange(v1.shape[1])])

def remez(func, interval, degree, error=None, maxiter=30, float_type=numpy.float128):
    """
        The remez algorithm is an iterative algorithm for finding the optimal polynomial for a giving function on a
    closed interval.
        Chebyshev showed that such a polynomial 'exists' and is 'unique', and meets the following:
            - If R(x) is a polynomial of degree N, then there are N+2 unknowns:
                the N+1 coefficients of the polynomial, and maximal value of the error function.
            - The error function has N+1 roots, and N+2 extrema (minima and maxima).
            - The extrema alternate in sign, and all have the same magnitude.
        The key to finding this polynomial is locating those locations withing then closed interval, that meets all
        three of these properties.
    If we know the location of the extrema of the error function, then we can write N+2 simultaneous equations:
        R(xi) + (-1)iE = f(xi)
    where E is the maximal error term, and xi are the abscissa values of the N+2 extrema of the error function.
    It is then trivial to solve the simultaneous equations to obtain the polynomial coefficients and the error term.
    Unfortunately we don't know where the extrema of the error function are located!

    The remez method is used to locate (hopefully converge in a timely manner) on such locations.

    1) Start by a 'good' estimate, using Chebyshev roots as the points in questions.
    note: this are only applicable on the interval [-1, 1], hence the Chebyshev roots need to be linearly mapped
        to the giving interval [a, b].
    2) Using polynomial interpolation or any other method to locate the initial set of coefficients ...
    3) Locate all local extrema there should N+2 such locations see: get_extrema
    4) create a new solution, (coefficients + error_term) using the extrema(s), if the error_term doesn't change
        by a certain amount quit since progress can no long er be made
        otherwise use the previous extrema(s) as the new locations and repeat steps 3, 4 ...
    """
    f = func if type(func) is numpy.ufunc else numpy.vectorize(func)  # vectorized non-numpy functions ...
    # numpy.pi is a float64 value, this should give us a bit more accuracy ...
    one, two, four, five, sixteen = imap(float_type, (1, 2, 4, 5, 16))
    pi = sixteen * numpy.arctan(one / five) - four * numpy.arctan(one / float_type(239))
    chebyshev_nodes = numpy.cos(  # locate all needed chebyshev nodes ...
        (((two * degree + one - two * numpy.arange(0, degree + 1, dtype=float_type)) * pi)/(two * degree + two))
    )
    # linearly map chebyshev nodes from (-1, 1) to the giving interval, scale + offset ...
    x = (numpy.diff(interval) / two) * chebyshev_nodes + numpy.sum(interval) / two
    fx = f(x)
    coefficients = solve(numpy.vander(x), fx)  # solve the system ...
    # relative error function .. bind the current coefficients to it ...
    rel_error_func = lambda v, coefficients=coefficients, f=f: (numpy.polyval(coefficients, v) - f(v))/f(v)
    alternating_sign = alternating_signs((degree + 2,))
    delta_error_term, error_term = 10, 1000
    x = remez_get_extremas(rel_error_func, interval, roots=x)  # get extremas from Chebyshev roots and use them for sol
    error = numpy.finfo(x.dtype).eps if error is None else error  # set the error to the floats machine epsilon ...
    while abs(delta_error_term) > error and maxiter:  # continue while making progress
        x = remez_get_extremas(
            lambda v, coefficients=coefficients, f=f: rel_error_func(v, coefficients, f), interval, x, accuracy=error
        )
        fx = f(x)
        new_solution = solve(  # solve the system of N + 2 equations to get a new solution and error term
            numpy.append(numpy.vander(x, degree + 1), (alternating_sign * numpy.abs(fx)).reshape(-1, 1), axis=1), fx
        )  # I think f(xi)*-1**i has to be added as the last term (E) in order for errorfunc to equioscillate at extrema
        delta_error_term = new_solution[-1] - error_term
        coefficients, error_term = new_solution[:-1], new_solution[-1]
        maxiter -= 1
    return coefficients

def test_score_mean_mt_mse():
    V3_true = np.vander(np.arange(3))
    V3_pred = np.vander(np.arange(1, 4))
    assert test.score_mean_mt_mse(V3_true, V3_pred) - 4.22222 < 1e-5

    V3_true_ma = np.ma.MaskedArray(V3_true)
    V3_true_ma.mask = np.zeros((3, 3))
    V3_true_ma.mask[2, :] = 1
    assert test.score_mean_mt_mse(V3_true_ma, V3_pred) == 2

 def test_basic(self):
     c = np.array([0, 1, -2, 3])
     v = vander(c)
     powers = np.array([[ 0,  0, 0,  0, 1],
                        [ 1,  1, 1,  1, 1],
                        [16, -8, 4, -2, 1],
                        [81, 27, 9,  3, 1]])
     # Check default value of N:
     yield (assert_array_equal, v, powers[:, 1:])
     # Check a range of N values, including 0 and 5 (greater than default)
     m = powers.shape[1]
     for n in range(6):
         v = vander(c, N=n)
         yield (assert_array_equal, v, powers[:, m-n:m])

def polyfit(x, y, deg, rcond=None, full=False):
    """%s

    Notes
    -----
        Any masked values in x is propagated in y, and vice-versa.
    """
    order = int(deg) + 1
    x = asarray(x)
    mx = getmask(x)
    y = asarray(y)
    if y.ndim == 1:
        m = mask_or(mx, getmask(y))
    elif y.ndim == 2:
        y = mask_rows(y)
        my = getmask(y)
        if my is not nomask:
            m = mask_or(mx, my[:,0])
        else:
            m = mx
    else:
        raise TypeError,"Expected a 1D or 2D array for y!"
    if m is not nomask:
        x[m] = y[m] = masked
    # Set rcond
    if rcond is None :
        if x.dtype in (np.single, np.csingle):
            rcond = len(x)*_single_eps
        else :
            rcond = len(x)*_double_eps
    # Scale x to improve condition number
    scale = abs(x).max()
    if scale != 0 :
        x = x / scale
    # solve least squares equation for powers of x
    v = vander(x, order)
    c, resids, rank, s = _lstsq(v, y.filled(0), rcond)
    # warn on rank reduction, which indicates an ill conditioned matrix
    if rank != order and not full:
        warnings.warn("Polyfit may be poorly conditioned", np.RankWarning)
    # scale returned coefficients
    if scale != 0 :
        if c.ndim == 1 :
            c /= np.vander([scale], order)[0]
        else :
            c /= np.vander([scale], order).T
    if full :
        return c, resids, rank, s, rcond
    else :
        return c

def test_graph_laplacian():
    for mat in (np.arange(10) * np.arange(10)[:, np.newaxis],
                np.ones((7, 7)),
                np.eye(19),
                np.vander(np.arange(4)) + np.vander(np.arange(4)).T,):
        sp_mat = sparse.csr_matrix(mat)
        for normed in (True, False):
            laplacian = graph_laplacian(mat, normed=normed)
            n_nodes = mat.shape[0]
            if not normed:
                np.testing.assert_array_almost_equal(laplacian.sum(axis=0),
                                                     np.zeros(n_nodes))
            np.testing.assert_array_almost_equal(laplacian.T, laplacian)
            np.testing.assert_array_almost_equal(
                laplacian, graph_laplacian(sp_mat, normed=normed).toarray())

def linest(*args, **kwargs):

    Y = args[0]
    X = args[1]
    
    if len(args) == 3:
        const = args[2]
        if isinstance(const,str):
            const = (const.lower() == "true")
    else:
        const = True
        
    degree = kwargs.get('degree',1)
    
    # build the vandermonde matrix
    A = np.vander(X, degree+1)
    
    if not const:
        # force the intercept to zero
        A[:,-1] = np.zeros((1,len(X)))
    
    # perform the fit
    (coefs, residuals, rank, sing_vals) = np.linalg.lstsq(A, Y)
        
    return coefs

    def test_dtypes(self):
        c = array([11, -12, 13], dtype=np.int8)
        v = vander(c)
        expected = np.array([[121,  11, 1],
                             [144, -12, 1],
                             [169,  13, 1]])
        yield (assert_array_equal, v, expected)

        c = array([1.0+1j, 1.0-1j])
        v = vander(c, N=3)
        expected = np.array([[ 2j, 1+1j, 1],
                             [-2j, 1-1j, 1]])
        # The data is floating point, but the values are small integers,
        # so assert_array_equal *should* be safe here (rather than, say,
        # assert_array_almost_equal).
        yield (assert_array_equal, v, expected)

def polybkdfit(q,sq,porder):
    '''
    Module to fit a polynomial background to s(q).
    
    Inputs: array q, array s(q), and desired order of the background polynomial to be fit.
    
    Returns: polynomial coefficients in array p.
    '''
    qscale=q[-1]
    qsc = q / qscale
    Mv0 = np.vander(qsc,porder+1)
    Mv1 = Mv0[:,:-1]
    yfq = q * (sq - 1.0)
    p,resids,rank,s=lstsq(Mv1,yfq)
    p /= np.vander([qscale],porder+1)[0,:-1]
    return p

def linest(*args, **kwargs): # Excel reference: https://support.office.com/en-us/article/LINEST-function-84d7d0d9-6e50-4101-977a-fa7abf772b6d

    Y = args[0].values()
    X = args[1].values()
    
    if len(args) == 3:
        const = args[2]
        if isinstance(const,str):
            const = (const.lower() == "true")
    else:
        const = True
        
    degree = kwargs.get('degree',1)
    
    # build the vandermonde matrix
    A = np.vander(X, degree+1)
    
    if not const:
        # force the intercept to zero
        A[:,-1] = np.zeros((1,len(X)))
    
    # perform the fit
    (coefs, residuals, rank, sing_vals) = np.linalg.lstsq(A, Y)
        
    return coefs

def fit_continuum(x, y, ivars, order=6, nsigma=[0.3,3.0], maxniter=50):
    """Fit the continuum using sigma clipping

    Args:
        x: The wavelengths
        y: The log-fluxes
        order: The polynomial order to use
        nsigma: The sigma clipping threshold: tuple (low, high)
        maxniter: The maximum number of iterations to do

    Returns:
        The value of the continuum at the wavelengths in x

    """
    A = np.vander(x - np.nanmean(x), order+1)
    m = np.ones(len(x), dtype=bool)
    for i in range(maxniter):
        m[ivars == 0] = 0  # mask out the bad pixels
        w = np.linalg.solve(np.dot(A[m].T, A[m]), np.dot(A[m].T, y[m]))
        mu = np.dot(A, w)
        resid = y - mu
        sigma = np.sqrt(np.nanmedian(resid**2))
        #m_new = np.abs(resid) < nsigma*sigma
        m_new = (resid > -nsigma[0]*sigma) & (resid < nsigma[1]*sigma)
        if m.sum() == m_new.sum():
            m = m_new
            break
        m = m_new
    return mu

    def __call__(self, xnew):
        saveshape = np.shape(xnew)
        xnew = np.ravel(xnew)
        res = np.empty_like(xnew)
        mask = (self.a <= xnew) & (xnew <= self.b)
        res[~mask] = self.fill
        xx = xnew.compress(mask)
        indxs = np.searchsorted(self.breaks[:-1], xx) - 1
        indxs = indxs.clip(0, len(self.breaks))
        pp = self.coeffs
        dx = xx - self.breaks.take(indxs)
        if True:
            v = pp[0, indxs]
            for i in xrange(1, self.order):
                v = dx * v + pp[i, indxs]
            values = v
        else:
            V = np.vander(dx, N=self.order)
            # values = np.diag(dot(V,pp[:,indxs]))
            dot = np.dot
            values = np.array([dot(V[k, :], pp[:, indxs[k]])
                              for k in xrange(len(xx))])

        res[mask] = values
        res.shape = saveshape
        return res

def polynomial_expansion_Vandermonde(a, degree=1):
    """
    Performs polynomial expansion of given :math:`a` array using *Vandermonde*
    method.

    Parameters
    ----------
    a : array_like
        :math:`a` array to expand.
    degree : int, optional
        Expanded polynomial degree.

    Returns
    -------
    ndarray
        Expanded :math:`a` array.

    References
    ----------
    :cite:`Wikipedia2003e`

    Examples
    --------
    >>> RGB = np.array([0.17224810, 0.09170660, 0.06416938])
    >>> polynomial_expansion_Vandermonde(RGB)  # doctest: +ELLIPSIS
    array([ 0.1722481 ,  0.0917066 ,  0.06416938,  1.        ])
    """

    a = as_float_array(a)

    a_e = np.transpose(np.vander(np.ravel(a), degree + 1))
    a_e = np.hstack(a_e.reshape(a_e.shape[0], -1, 3))

    return np.squeeze(a_e[:, 0:a_e.shape[-1] - a.shape[-1] + 1])

def lomb_scargle_estimator(x, y, yerr=None,
                           min_period=None, max_period=None,
                           filter_period=None,
                           max_peaks=2,
                           **kwargs):
    """Estimate period of a time series using the periodogram

    Args:
        x (ndarray[N]): The times of the observations
        y (ndarray[N]): The observations at times ``x``
        yerr (Optional[ndarray[N]]): The uncertainties on ``y``
        min_period (Optional[float]): The minimum period to consider
        max_period (Optional[float]): The maximum period to consider
        filter_period (Optional[float]): If given, use a high-pass filter to
            down-weight period longer than this
        max_peaks (Optional[int]): The maximum number of peaks to return
            (default: 2)

    Returns:
        A dictionary with the computed ``periodogram`` and the parameters for
        up to ``max_peaks`` peaks in the periodogram.

    """
    if min_period is not None:
        kwargs["maximum_frequency"] = 1.0 / min_period
    if max_period is not None:
        kwargs["minimum_frequency"] = 1.0 / max_period

    # Estimate the power spectrum
    model = LombScargle(x, y, yerr)
    freq, power = model.autopower(method="fast", normalization="psd", **kwargs)
    power /= len(x)
    power_est = np.array(power)

    # Filter long periods
    if filter_period is not None:
        freq0 = 1.0 / filter_period
        filt = 1.0 / np.sqrt(1 + (freq0 / freq) ** (2*3))
        power *= filt

    # Find and fit peaks
    peak_inds = (power[1:-1] > power[:-2]) & (power[1:-1] > power[2:])
    peak_inds = np.arange(1, len(power)-1)[peak_inds]
    peak_inds = peak_inds[np.argsort(power[peak_inds])][::-1]
    peaks = []
    for i in peak_inds[:max_peaks]:
        A = np.vander(freq[i-1:i+2], 3)
        w = np.linalg.solve(A, np.log(power[i-1:i+2]))
        sigma2 = -0.5 / w[0]
        freq0 = w[1] * sigma2
        peaks.append(dict(
            log_power=w[2] + 0.5*freq0**2 / sigma2,
            period=1.0 / freq0,
            period_uncert=np.sqrt(sigma2 / freq0**4),
        ))

    return dict(
        periodogram=(freq, power_est),
        peaks=peaks,
    )

def get_Ab(n,m):
	N = n + 1 				# N = n + 1
	X = np.linspace(-5,5,m)
	A = np.vander(X,N) # Sehr zu empfehlender Trick

	F = lambda x: 1 / (1 + x**2)
	return A, F(X)