Python scipy.linalg.toeplitz() Examples

def gaussian_emd(x, y, sigma=1.0, distance_scaling=1.0):
    ''' Gaussian kernel with squared distance in exponential term replaced by EMD
    Args:
      x, y: 1D pmf of two distributions with the same support
      sigma: standard deviation
    '''
    support_size = max(len(x), len(y))
    d_mat = toeplitz(range(support_size)).astype(np.float)
    distance_mat = d_mat / distance_scaling

    # convert histogram values x and y to float, and make them equal len
    x = x.astype(np.float)
    y = y.astype(np.float)
    if len(x) < len(y):
        x = np.hstack((x, [0.0] * (support_size - len(x))))
    elif len(y) < len(x):
        y = np.hstack((y, [0.0] * (support_size - len(y))))

    emd = pyemd.emd(x, y, distance_mat)
    return np.exp(-emd * emd / (2 * sigma * sigma)) 

def toepz(self):
        cormat = np.zeros((self.nkdim, self.nkdim))

        epsilon = (1 - np.max(self.rho)) / (1 + np.max(self.rho)) - .01

        for i in np.arange(self.k):
            t = np.insert(np.power(self.rho[i], np.arange(1, self.nk[i])), 0, 1)
            cor = toeplitz(t)
            if i == 0:
                cormat[0:self.nk[0], 0:self.nk[0]] = cor
            if i != 0:
                cormat[np.sum(self.nk[0:i]):np.sum(self.nk[0:i + 1]),
                np.sum(self.nk[0:i]):np.sum(self.nk[0:i + 1])] = cor

        np.fill_diagonal(cormat, 1 - epsilon)

        cormat = self._generate_noise(cormat, self.nkdim, self.M, epsilon)

        return cormat 

def hub(self):
        cormat = np.zeros((self.nkdim, self.nkdim))

        for i in np.arange(self.k):
            cor = toeplitz(self._fill_hub_matrix(self.rho[i, 0], self.rho[i, 1], self.power, self.nk[i]))
            if i == 0:
                cormat[0:self.nk[0], 0:self.nk[0]] = cor
            if i != 0:
                cormat[np.sum(self.nk[0:i]):np.sum(self.nk[0:i + 1]),
                np.sum(self.nk[0:i]):np.sum(self.nk[0:i + 1])] = cor
            tau = (np.max(self.rho[i]) - np.min(self.rho[i])) / (self.nk[i] - 2)

        epsilon = 0.08 #(1 - np.min(rho) - 0.75 * np.min(tau)) - 0.01

        np.fill_diagonal(cormat, 1 - epsilon)

        cormat = self._generate_noise(cormat, self.nkdim, self.M, epsilon)

        return cormat 

def corr_ar(k_vars, ar):
    '''create autoregressive correlation matrix

    This might be MA, not AR, process if used for residual process - check

    Parameters
    ----------
    ar : array_like, 1d
        AR lag-polynomial including 1 for lag 0


    '''
    from scipy.linalg import toeplitz
    if len(ar) < k_vars:
        ar_ = np.zeros(k_vars)
        ar_[:len(ar)] = ar
        ar = ar_

    return toeplitz(ar) 

def toeplitz_mul_v(toep,v):
    '''multiply a toeplitz matrix A by vector v.

    Args:
        toep (ndarray): representation fo multilevel toeplitz matrix, i.e.
            the first column of A in proper shape.
        v (ndarray): vector to be multiplied. Should be reshaped to the same shape
             as toep. Should be the same reshape order (column first/row first) as circ.
    Returns:
        result of multiplication.
    '''
    
    #xp = cupy.get_array_module(toep)
    circ,tmpv = embed_toep2circ(toep,v)
    res = circ_mul_v(circ,tmpv)
    slices = [slice(None)]*res.ndim
    slices[-1] = slice(0,v.shape[-1])
    slices[-2] = slice(0,v.shape[-2])
    return(res[tuple(slices)]) 

def corr_arma(k_vars, ar, ma):
    '''create arma correlation matrix

    converts arma to autoregressive lag-polynomial with k_var lags

    ar and arma might need to be switched for generating residual process

    Parameters
    ----------
    ar : array_like, 1d
        AR lag-polynomial including 1 for lag 0
    ma : array_like, 1d
        MA lag-polynomial

    '''
    from scipy.linalg import toeplitz
    from statsmodels.tsa.arima_process import arma2ar

    ar = arma2ar(ar, ma, nobs=k_vars)[:k_vars]  #bug in arma2ar

    return toeplitz(ar) 

def _params2cov(self, params, nobs):
        '''get autocovariance matrix from ARMA regression parameter

        ar parameters are assumed to have rhs parameterization

        '''
        ar = np.r_[[1], -params[:self.nar]]
        ma = np.r_[[1], params[-self.nma:]]
        #print('ar', ar
        #print('ma', ma
        #print('nobs', nobs
        autocov = arma_acovf(ar, ma, nobs=nobs)
        #print('arma_acovf(%r, %r, nobs=%d)' % (ar, ma, nobs)
        #print(autocov.shape
        #something is strange  fixed in aram_acovf
        autocov = autocov[:nobs]
        sigma = toeplitz(autocov)
        return sigma 

def setup_class(cls):
        from .results.results_regression import LongleyGls

        data = longley.load()
        exog = add_constant(np.column_stack(
            (data.exog[:, 1], data.exog[:, 4])), prepend=False)
        tmp_results = OLS(data.endog, exog).fit()
        rho = np.corrcoef(tmp_results.resid[1:],
                          tmp_results.resid[:-1])[0][1]  # by assumption
        order = toeplitz(np.arange(16))
        sigma = rho**order
        GLS_results = GLS(data.endog, exog, sigma=sigma).fit()
        cls.res1 = GLS_results
        cls.res2 = LongleyGls()
        # attach for test_missing
        cls.sigma = sigma
        cls.exog = exog
        cls.endog = data.endog 

def gaussian_emd(x, y, sigma=1.0, distance_scaling=1.0):
    ''' Gaussian kernel with squared distance in exponential term replaced by EMD
    Args:
      x, y: 1D pmf of two distributions with the same support
      sigma: standard deviation
    '''
    support_size = max(len(x), len(y))
    d_mat = toeplitz(range(support_size)).astype(np.float)
    distance_mat = d_mat / distance_scaling

    # convert histogram values x and y to float, and make them equal len
    x = x.astype(np.float)
    y = y.astype(np.float)
    if len(x) < len(y):
        x = np.hstack((x, [0.0] * (support_size - len(x))))
    elif len(y) < len(x):
        y = np.hstack((y, [0.0] * (support_size - len(y))))

    emd = pyemd.emd(x, y, distance_mat)
    return np.exp(-emd * emd / (2 * sigma * sigma)) 

def Rmtx_ri(coef_ri, K, D, L):
    ''' Split T matrix in rea/imaginary representation '''
    coef_ri = np.squeeze(coef_ri)
    coef_r = coef_ri[:K + 1]
    coef_i = coef_ri[K + 1:]
    R_r = linalg.toeplitz(np.concatenate((np.array([coef_r[-1]]),
                                          np.zeros(L - K - 1))),
                          np.concatenate((coef_r[::-1],
                                          np.zeros(L - K - 1)))
                          )
    R_i = linalg.toeplitz(np.concatenate((np.array([coef_i[-1]]),
                                          np.zeros(L - K - 1))),
                          np.concatenate((coef_i[::-1],
                                          np.zeros(L - K - 1)))
                          )
    return np.dot(np.vstack((np.hstack((R_r, -R_i)), np.hstack((R_i, R_r)))), D) 

def corr_ar(k_vars, ar):
    '''create autoregressive correlation matrix

    This might be MA, not AR, process if used for residual process - check

    Parameters
    ----------
    ar : array_like, 1d
        AR lag-polynomial including 1 for lag 0


    '''
    from scipy.linalg import toeplitz
    if len(ar) < k_vars:
        ar_ = np.zeros(k_vars)
        ar_[:len(ar)] = ar
        ar = ar_

    return toeplitz(ar) 

def Tmtx_ri(b_ri, K, D, L):
    """
    build convolution matrix associated with b_ri

    Parameters
    ----------
    b_ri: 
        a real-valued vector
    K: 
        number of Diracs
    D1: 
        expansion matrix for the real-part
    D2: 
        expansion matrix for the imaginary-part
    """
    b_ri = np.dot(D, b_ri)
    b_r = b_ri[:L]
    b_i = b_ri[L:]
    Tb_r = linalg.toeplitz(b_r[K:], b_r[K::-1])
    Tb_i = linalg.toeplitz(b_i[K:], b_i[K::-1])
    return np.vstack((np.hstack((Tb_r, -Tb_i)), np.hstack((Tb_i, Tb_r)))) 

def test_solve_equivalence():
    # For toeplitz matrices, solve_toeplitz() should be equivalent to solve().
    random = np.random.RandomState(1234)
    for n in (1, 2, 3, 10):
        c = random.randn(n)
        if random.rand() < 0.5:
            c = c + 1j * random.randn(n)
        r = random.randn(n)
        if random.rand() < 0.5:
            r = r + 1j * random.randn(n)
        y = random.randn(n)
        if random.rand() < 0.5:
            y = y + 1j * random.randn(n)

        # Check equivalence when both the column and row are provided.
        actual = solve_toeplitz((c,r), y)
        desired = solve(toeplitz(c, r=r), y)
        assert_allclose(actual, desired)

        # Check equivalence when the column is provided but not the row.
        actual = solve_toeplitz(c, b=y)
        desired = solve(toeplitz(c), y)
        assert_allclose(actual, desired) 

def test_unstable():
    # this is a "Gaussian Toeplitz matrix", as mentioned in Example 2 of
    # I. Gohbert, T. Kailath and V. Olshevsky "Fast Gaussian Elimination with
    # Partial Pivoting for Matrices with Displacement Structure"
    # Mathematics of Computation, 64, 212 (1995), pp 1557-1576
    # which can be unstable for levinson recursion.

    # other fast toeplitz solvers such as GKO or Burg should be better.
    random = np.random.RandomState(1234)
    n = 100
    c = 0.9 ** (np.arange(n)**2)
    y = random.randn(n)

    solution1 = solve_toeplitz(c, b=y)
    solution2 = solve(toeplitz(c), y)

    assert_allclose(solution1, solution2) 

def non_toeplitz_covariance(data, window_size):
    """
    Get scaled non- Toeplitz covariance matrix, which may be able to account
    for non-stationary data-errors.

    Parameters
    ----------
    data : :class:`numpy.ndarray`
        of data to estimate non-stationary error matrix for
    window_size : int
        samples to take on running rmse estimation over data

    Returns
    -------
    :class:`numpy.ndarray` (size data, size data)
    """
    toeplitz, stds = toeplitz_covariance(data, window_size)
    return toeplitz * stds[:, num.newaxis] * stds[num.newaxis, :] 

def corr_arma(k_vars, ar, ma):
    '''create arma correlation matrix

    converts arma to autoregressive lag-polynomial with k_var lags

    ar and arma might need to be switched for generating residual process

    Parameters
    ----------
    ar : array_like, 1d
        AR lag-polynomial including 1 for lag 0
    ma : array_like, 1d
        MA lag-polynomial

    '''
    from scipy.linalg import toeplitz
    from statsmodels.tsa.arima_process import arma2ar

    ar = arma2ar(ar, ma, nobs=k_vars)[:k_vars]  #bug in arma2ar

    return toeplitz(ar) 

def estimate(self, nbcorr=np.nan, numpsd=-1):
        fft_length, _ = self.check_params()
        if np.isnan((nbcorr)):
            nbcorr = self.ordar

        # -------- estimate correlation from psd
        full_psd = self.psd[numpsd]
        full_psd = np.c_[full_psd, np.conjugate(full_psd[:, :0:-1])]
        correl = fftpack.ifft(full_psd[0], fft_length, 0).real

        # -------- estimate AR part
        col1 = correl[self.ordma:self.ordma + nbcorr]
        row1 = correl[np.abs(
            np.arange(self.ordma, self.ordma - self.ordar, -1))]
        R = linalg.toeplitz(col1, row1)
        r = -correl[self.ordma + 1:self.ordma + nbcorr + 1]
        AR = linalg.solve(R, r)
        self.AR_ = AR

        # -------- estimate correlation of MA part

        # -------- estimate MA part
        if self.ordma == 0:
            sigma2 = correl[0] + np.dot(AR, correl[1:self.ordar + 1])
            self.MA = np.ones(1) * np.sqrt(sigma2)
        else:
            raise NotImplementedError(
                'arma: estimation of the MA part not yet implemented') 

def Rmtx_ri(coef_ri, K, D, L):
    coef_ri = np.squeeze(coef_ri)
    coef_r = coef_ri[:K + 1]
    coef_i = coef_ri[K + 1:]
    R_r = linalg.toeplitz(np.concatenate((np.array([coef_r[-1]]),
                                          np.zeros(L - K - 1))),
                          np.concatenate((coef_r[::-1],
                                          np.zeros(L - K - 1)))
                          )
    R_i = linalg.toeplitz(np.concatenate((np.array([coef_i[-1]]),
                                          np.zeros(L - K - 1))),
                          np.concatenate((coef_i[::-1],
                                          np.zeros(L - K - 1)))
                          )
    return np.dot(np.vstack((np.hstack((R_r, -R_i)), np.hstack((R_i, R_r)))), D) 

def yule_walker_acov(acov, order=1, method="unbiased", df=None, inv=False):
    """
    Estimate AR(p) parameters from acovf using Yule-Walker equation.


    Parameters
    ----------
    acov : array-like, 1d
        auto-covariance
    order : integer, optional
        The order of the autoregressive process.  Default is 1.
    inv : bool
        If inv is True the inverse of R is also returned.  Default is False.

    Returns
    -------
    rho : ndarray
        The estimated autoregressive coefficients
    sigma
        TODO
    Rinv : ndarray
        inverse of the Toepliz matrix

    """


    R = toeplitz(r[:-1])

    rho = np.linalg.solve(R, r[1:])
    sigmasq = r[0] - (r[1:]*rho).sum()
    if inv == True:
        return rho, np.sqrt(sigmasq), np.linalg.inv(R)
    else:
        return rho, np.sqrt(sigmasq) 

def get_ar_coef(y, sn, p, add_lag, pad=None):
    if add_lag is 'p':
        max_lag = p * 2
    else:
        max_lag = p + add_lag
    cov = acovf(y, fft=True)
    C_mat = toeplitz(cov[:max_lag], cov[:p]) - sn**2 * np.eye(max_lag, p)
    g = lstsq(C_mat, cov[1:max_lag + 1])[0]
    if pad:
        res = np.zeros(pad)
        res[:len(g)] = g
        return res
    else:
        return g 

def coef_restriction_diffseq(n_coeffs, degree=1, n_vars=None, position=0, base_idx=0):
    #check boundaries, returns "valid" ?

    if degree == 1:
        diff_coeffs = [-1, 1]
        n_points = 2
    elif degree > 1:
        from scipy import misc
        n_points = next_odd(degree + 1)  #next odd integer after degree+1
        diff_coeffs = misc.central_diff_weights(n_points, ndiv=degree)

    dff = np.concatenate((diff_coeffs, np.zeros(n_coeffs - len(diff_coeffs))))
    from scipy import linalg
    reduced = linalg.toeplitz(dff, np.zeros(n_coeffs - len(diff_coeffs) + 1)).T
    #reduced = np.kron(np.eye(n_coeffs-n_points), diff_coeffs)

    if n_vars is None:
        return reduced
    else:
        full = np.zeros((n_coeffs-1, n_vars))
        full[:, position:position+n_coeffs] = reduced
        return full


##
##    R = np.c_[np.zeros((n_groups, k_vars-1)), np.eye(n_groups)]
##    r = np.zeros(n_groups)
##    R = np.c_[np.zeros((n_groups-1, k_vars)),
##              np.eye(n_groups-1)-1./n_groups * np.ones((n_groups-1, n_groups-1))] 

def covariance_matrix_grid(self, endog_expval, index):

        from scipy.linalg import toeplitz
        r = np.zeros(len(endog_expval))
        r[0] = 1
        r[1:self.max_lag + 1] = self.dep_params
        return toeplitz(r), True 

def embed_toep2circ(toep,v=None):
    '''embed toeplitz matrix to circulant matrix.

    Args:
        toep (ndarray): representation of multilevel toeplitz matrix, i.e.
            the first column of A in proper shape.
        v (ndarray): embed a vector together with toep.
    Returns:
        representation of embedded multilevel circulant matrix and embedded vector.
    '''
    if use_gpu > 0:
        import cupy
        xp = cupy.get_array_module(toep)
    else:
        xp = np
        
#    circ = xp.zeros((2*xp.array(toep.shape)).astype(xp.int))
    circ = xp.zeros((2*toep.shape[0],2*toep.shape[1]))
    s = []
    for idim in range(toep.ndim):
        s.append(slice(0,toep.shape[idim]))
    circ[tuple(s)] = toep
    if not v is None:
        if v.ndim == toep.ndim:
            resv = xp.zeros_like(circ)
            resv[tuple(s)] = v
        else:
            resv = xp.zeros((v.shape[0],*circ.shape))
            resv[tuple(slice(None),*s)] = v
    for idim in range(toep.ndim):
        s[idim] = slice(toep.shape[idim]+1,None)
        s2 = s.copy()
        s2[idim] = slice(toep.shape[idim]-1,0,-1)
        circ[tuple(s)] = circ[tuple(s2)]
        s[idim] = slice(None)
    if v is None:
        return circ
    else:
        return circ,resv 

def _compute_reference_correlations(reference_sources, filters_len):
    """Compute the inner products between delayed versions of reference_sources
    reference is nsrc X nsamp X nchan.
    Returns
    * G, matrix : nsrc X nsrc X nchan X nchan X filters_len X filters_len
    * sf, reference spectra: nsrc X nchan X filters_len"""

    # reshape references as nsrc X nchan X nsampl
    (nsrc, nsampl, nchan) = reference_sources.shape
    reference_sources = np.moveaxis(reference_sources, (1), (2))

    # zero padding and FFT of references
    reference_sources = _zeropad(reference_sources, filters_len - 1, axis=2)
    n_fft = int(2**np.ceil(np.log2(nsampl + filters_len - 1.)))
    sf = scipy.fftpack.fft(reference_sources, n=n_fft, axis=2)

    # compute intercorrelation between sources
    G = np.zeros((nsrc, nsrc, nchan, nchan, filters_len, filters_len))
    for ((i, c1), (j, c2)) in itertools.combinations_with_replacement(
        itertools.product(
            list(range(nsrc)), list(range(nchan))
        ),
        2
    ):

        ssf = sf[j, c2] * np.conj(sf[i, c1])
        ssf = np.real(scipy.fftpack.ifft(ssf))
        ss = toeplitz(
            np.hstack((ssf[0], ssf[-1:-filters_len:-1])),
            r=ssf[:filters_len]
        )
        G[j, i, c2, c1] = ss
        G[i, j, c1, c2] = ss.T
    return G, sf 

def __getitem__(self, key):
        slc0, slc1 = self._unpack_index(key)
        i0, i1 = self._process_slice(slc0, self.shape[0])
        j0, j1 = self._process_slice(slc1, self.shape[1])
        C, R = self._c, self._r

        # symmetric query
        if (i0 == j0) and (i1 == j1):
            c = C[0 : (i1 - i0)]
            r = R[0 : (j1 - j0)]

        # asymmetric query
        else:
            transpose = False
            # tril
            if j0 < i0 or (i0 == j0 and i1 < j1):
                # tranpose the matrix, query,
                # then transpose the result
                i0, i1, j0, j1 = j0, j1, i0, i1
                C, R = R, C
                transpose = True

            c = np.r_[R[(j0 - i0) : max(0, j0 - i1) : -1], C[0 : max(0, i1 - j0)]]
            r = R[(j0 - i0) : (j1 - i0)]

            if transpose:
                c, r = r, c

        return toeplitz(c, r) 

def make_cis_obsexp_fetcher(clr, expected, weight_name="weight"):
    """
    Construct a function that returns intra-chromosomal OBS/EXP.

    Parameters
    ----------
    clr : cooler.Cooler
        Observed matrix.
    expected : tuple of (DataFrame, str)
        Diagonal summary statistics for each chromosome, and name of the column
        to use
    weight_name : str
        Name of the column in the clr.bins to use as balancing weights

    Returns
    -------
    getexpected(chrom, _). 2nd arg is ignored.

    """
    expected, name = expected
    expected = {k: x.values for k, x in expected.groupby("chrom")[name]}

    def _fetch_cis_oe(reg1, reg2):
        obs_mat = clr.matrix(balance=weight_name).fetch(reg1)
        exp_mat = toeplitz(expected[reg1[0]][: obs_mat.shape[0]])
        return obs_mat / exp_mat

    return _fetch_cis_oe 

def timeMatrix(n):
	col = [0] * n
	col[0] = 1
	raw = [0] * n
	raw[0] = 1
	raw[1] = -1
	return matrix(toeplitz(col, raw)) 

def timeMatrix(n):
	col = [0] * n
	col[0] = 1
	raw = [0] * n
	raw[0] = 1
	raw[1] = -1
	return matrix(toeplitz(col, raw)) 

def estimate(self, nbcorr=np.nan, numpsd=-1):
        fft_length, _ = self.check_params()
        if np.isnan((nbcorr)):
            nbcorr = self.ordar

        # -------- estimate correlation from psd
        full_psd = self.psd[numpsd]
        full_psd = np.c_[full_psd, np.conjugate(full_psd[:, :0:-1])]
        correl = fftpack.ifft(full_psd[0], fft_length, 0).real

        # -------- estimate AR part
        col1 = correl[self.ordma:self.ordma + nbcorr]
        row1 = correl[np.abs(
            np.arange(self.ordma, self.ordma - self.ordar, -1))]
        R = linalg.toeplitz(col1, row1)
        r = -correl[self.ordma + 1:self.ordma + nbcorr + 1]
        AR = linalg.solve(R, r)
        self.AR_ = AR

        # -------- estimate correlation of MA part

        # -------- estimate MA part
        if self.ordma == 0:
            sigma2 = correl[0] + np.dot(AR, correl[1:self.ordar + 1])
            self.MA = np.ones(1) * np.sqrt(sigma2)
        else:
            raise NotImplementedError(
                'arma: estimation of the MA part not yet implemented') 

def block_toep2_sym(a):
    '''generate full representation of 2-level symmetric toeplitz matrix

    Args:
        a (ndarray): 1-st column of the symmetrec toeplitz matrix in proper shape.
    Returns:
        Full filled toeplitz matrix.
    '''
    if use_gpu > 0:
        import cupy
        xp = cupy.get_array_module(a)
        if xp is cupy:
            a = xp.asnumpy(a)
    else:
        xp = np
        a = np.asnumpy(a)
        
    A1 = []
    n0,n1 = a.shape
    for i in range(n1):
        A1.append(splin.toeplitz(a[:,i]))
    A = np.empty((n1,n0,n1,n0))
    for i in range(n1):
        for j in range(n1):
            A[i,:,j,:] = A1[np.int(np.abs(i-j))]
    A.shape = (n0*n1,n0*n1)
    A = xp.asarray(A)
    
    return(A)