Python linalg.circulant方法代码示例

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def _process_masks(self, mask_obs, self_mask=False):
        '''
            mask_obs will be a (n_agent, n_object) boolean matrix. If the mask is over non-agent
                objects then we do nothing. If the mask is over other agents (self_mask is True),
                then we permute each row such that the mask is consistent with the circulant
                permutation used for self observations
        '''
        new_mask = mask_obs.copy()
        if self_mask:
            assert np.all(new_mask.shape == np.array((self.n_agents, self.n_agents)))
            # Permute each row to the right by one more than the previous
            # E.g., [[1,2],[3,4]] -> [[1,2],[4,3]]
            idx = circulant(np.arange(self.n_agents))
            new_mask = new_mask[np.arange(self.n_agents)[:, None], idx]
            new_mask = new_mask[:, 1:]  # Remove self observation
        return new_mask 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def toeplitz_strang_circ_approx(r, matrix=False):
    '''
    Circulant approximation to a symetric Toeplitz matrix
    by Gil Strang

    Parameters
    ----------
    r: ndarray
        the first row of the symmetric Toeplitz matrix
    matrix: bool, optional
        if True, the full symetric Toeplitz matrix is returned,
        otherwise, only the first column
    '''

    n = r.shape[0]
    c = r.copy()
    m = n // 2 if n % 2 == 0 else (n - 1) // 2
    c[-m:] = r[m:0:-1]

    if matrix:
        return la.circulant(c)
    else:
        return c 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def circ_mul_v(circ,v,eigs=None):
    ''' multiply a circulant matrix A by multi vector v.

    Args:
        circ (ndarray): representation of the multilevel circulant matrix A, i.e.
             the first column of A in proper shape.
        v (ndarray): vector to be multiplied. Should be reshaped to the same shape
             as circ. Should be the same reshape order (column first/row first) as circ.
    Returns:
        result of multiplication.
    '''
    if use_gpu > 0:
        import cupy
        xp = cupy.get_array_module(circ)
    else:
        xp = np
    
    if eigs is None:
        eigs = circ_eigs(circ)
    tmp = xp.real(xp.fft.ifft2(xp.fft.fft2(v,norm='ortho')*eigs,norm='ortho'))
    if xp is cupy:
        return tmp.astype(xp.float32)
    else:
        return tmp 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def observation(self, obs):
        new_obs = {}
        for k, v in obs.items():
            if 'mask' in k:
                new_obs[k] = self._process_masks(obs[k], self_mask=(k in self.keys_self))
            elif k in self.keys_self:
                new_obs[k + '_self'] = obs[k]
                new_obs[k] = obs[k][circulant(np.arange(self.n_agents))]
                new_obs[k] = new_obs[k][:, 1:, :]  # Remove self observation
            elif k in self.keys_copy:
                new_obs[k] = obs[k]
            else:
                new_obs[k] = np.tile(v, self.n_agents).reshape([v.shape[0], self.n_agents, v.shape[1]]).transpose((1, 0, 2))

        return new_obs 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def make_pd(start, n):
  """Deterministically create a positive definite matrix."""
  x = np.tril(linalg.circulant(np.arange(start, start + n)))
  return np.dot(x, x.T) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def getFDMatrix(N, order, dx):
  if order == 1:
    stencil = [-1.0, 1.0]
    zero_pos = 2
    coeff = 1.0

  elif order == 2:
    stencil = [1.0, -4.0, 3.0]
    coeff = 1.0 / 2.0
    zero_pos = 3

  elif order == 3:
    stencil = [1.0, -6.0, 3.0, 2.0]
    coeff = 1.0 / 6.0
    zero_pos = 3

  elif order == 4:
    stencil = [-5.0, 30.0, -90.0, 50.0, 15.0]
    coeff = 1.0 / 60.0
    zero_pos = 4

  elif order == 5:
    stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0]
    coeff = 1.0 / 60.0
    zero_pos = 5
  else:
    print("Order " + order + " not implemented.")

  first_col = np.zeros(N)

  # Because we need to specific first column (not row) in circulant, flip stencil array
  first_col[0:np.size(stencil)] = np.flipud(stencil)

  # Circulant shift of coefficient column so that entry number zero_pos becomes first entry
  first_col = np.roll(first_col, -np.size(stencil) + zero_pos, axis=0)

  return sp.csc_matrix(coeff * (1.0 / dx) * LA.circulant(first_col)) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def getUpwindMatrix(N, dx):
     
  #stencil    = [-1.0, 1.0]
  #zero_pos = 2
  #coeff      = 1.0
  
  #stencil    = [1.0, -4.0, 3.0]
  #coeff      = 1.0/2.0
  #zero_pos   = 3
  
  #stencil    = [1.0, -6.0, 3.0, 2.0]
  #coeff      = 1.0/6.0
  #zero_pos   = 3
  
  #stencil  = [-5.0, 30.0, -90.0, 50.0, 15.0]
  #coeff    = 1.0/60.0
  #zero_pos = 4
  
  stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0]
  coeff   = 1.0/60.0
  zero_pos = 5
  
  first_col = np.zeros(N)
  
  # Because we need to specific first column (not row) in circulant, flip stencil array
  first_col[0:np.size(stencil)] = np.flipud(stencil)

  # Circulant shift of coefficient column so that entry number zero_pos becomes first entry
  first_col = np.roll(first_col, -np.size(stencil)+zero_pos, axis=0)

  return sp.csc_matrix( coeff*(1.0/dx)*la.circulant(first_col) ) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def getHorizontalDx(N, dx, order):
    if order == 1:
        stencil = [-1.0, 1.0]
        zero_pos = 2
        coeff = 1.0

    elif order == 2:
        stencil = [1.0, -4.0, 3.0]
        coeff = 1.0 / 2.0
        zero_pos = 3

    elif order == 3:
        stencil = [1.0, -6.0, 3.0, 2.0]
        coeff = 1.0 / 6.0
        zero_pos = 3

    elif order == 4:
        stencil = [-5.0, 30.0, -90.0, 50.0, 15.0]
        coeff = 1.0 / 60.0
        zero_pos = 4

    elif order == 5:
        stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0]
        coeff = 1.0 / 60.0
        zero_pos = 5
    else:
        print("Order " + order + " not implemented.")

    first_col = np.zeros(N)

    # Because we need to specific first column (not row) in circulant, flip stencil array
    first_col[0:np.size(stencil)] = np.flipud(stencil)

    # Circulant shift of coefficient column so that entry number zero_pos becomes first entry
    first_col = np.roll(first_col, -np.size(stencil) + zero_pos, axis=0)

    return sp.csc_matrix(coeff * (1.0 / dx) * la.circulant(first_col)) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def getUpwindMatrix(N, dx, order):
    if order == 1:
        stencil = [-1.0, 1.0]
        zero_pos = 2
        coeff = 1.0

    elif order == 2:
        stencil = [1.0, -4.0, 3.0]
        coeff = 1.0 / 2.0
        zero_pos = 3

    elif order == 3:
        stencil = [1.0, -6.0, 3.0, 2.0]
        coeff = 1.0 / 6.0
        zero_pos = 3

    elif order == 4:
        stencil = [-5.0, 30.0, -90.0, 50.0, 15.0]
        coeff = 1.0 / 60.0
        zero_pos = 4

    elif order == 5:
        stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0]
        coeff = 1.0 / 60.0
        zero_pos = 5
    else:
        sys.exit('Order ' + str(order) + ' not implemented')

    first_col = np.zeros(N)

    # Because we need to specific first column (not row) in circulant, flip stencil array
    first_col[0:np.size(stencil)] = np.flipud(stencil)

    # Circulant shift of coefficient column so that entry number zero_pos becomes first entry
    first_col = np.roll(first_col, -np.size(stencil) + zero_pos, axis=0)

    return sp.csc_matrix(coeff * (1.0 / dx) * la.circulant(first_col)) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def _setup(self, config):
        torch.manual_seed(config['seed'])
        self.model = ButterflyProduct(size=config['size'],
                                      complex=False,
                                      fixed_order=config['fixed_order'],
                                      softmax_fn=config['softmax_fn'])
        if (not config['fixed_order']) and config['softmax_fn'] == 'softmax':
            self.semantic_loss_weight = config['semantic_loss_weight']
        self.optimizer = optim.Adam(self.model.parameters(), lr=config['lr'])
        self.n_steps_per_epoch = config['n_steps_per_epoch']
        size = config['size']
        # Need to transpose as dct acts on rows of matrix np.eye, not columns
        n = size
        np.random.seed(0)
        x = np.random.randn(n)
        C = circulant(x)
        self.target_matrix = torch.tensor(C, dtype=torch.float)
        arange_ = np.arange(size)
        dct_perm = np.concatenate((arange_[::2], arange_[::-2]))
        br_perm = bitreversal_permutation(size)
        assert config['perm'] in ['id', 'br', 'dct']
        if config['perm'] == 'id':
            self.perm = torch.arange(size)
        elif config['perm'] == 'br':
            self.perm = br_perm
        elif config['perm'] == 'dct':
            self.perm = torch.arange(size)[dct_perm][br_perm]
        else:
            assert False, 'Wrong perm in config' 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def toeplitz_opt_circ_approx(r, matrix=False):
    ''' 
    Optimal circulant approximation of a symmetric Toeplitz matrix
    by Tony F. Chan

    Parameters
    ----------
    r: ndarray
        the first row of the symmetric Toeplitz matrix
    matrix: bool, optional
        if True, the full symetric Toeplitz matrix is returned,
        otherwise, only the first column
    '''

    n = r.shape[0]

    r_rev = np.zeros(r.shape)
    r_rev[1:] = r[:0:-1]

    i = np.arange(n)
    c = (i * r_rev + (n - i) * r) / n
    c[1:] = c[:0:-1]

    if matrix:
        return la.circulant(c)
    else:
        return c 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def circ_eigs(circ,dtype=np.complex64):
    ''' calculate eigenvalues of multilevel circulant matrix A

    Args:
        circ (ndarray): representation of the multilevel circulant matrix A, i.e.
             the first column of A in proper shape.
    Returns:
        eigenvalues of A.
    '''
    if use_gpu > 0:
        import cupy
        xp = cupy.get_array_module(circ)
    else:
        xp = np
    return xp.fft.fft2(circ,norm='ortho').astype(dtype)*xp.sqrt(np.prod(circ.shape)) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
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 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def gradient(f, x=None, dx=1, axis=-1):
    """
    Return the gradient of 1 or 2-dimensional array.
    The gradient is computed using central differences in the interior
    and first differences at the boundaries.

    Irregular sampling is supported (it isn't supported by np.gradient)

    Parameters
    ----------
    f : 1d or 2d numpy array
        Input array.
    x : array_like, optional
       Points where the function f is evaluated. It must be of the same
       length as ``f.shape[axis]``.
       If None, regular sampling is assumed (see dx)
    dx : float, optional
       If `x` is None, spacing given by `dx` is assumed. Default is 1.
    axis : int, optional
       The axis along which the difference is taken.

    Returns
    -------
    out : array_like
        Returns the gradient along the given axis.

    Notes
    -----
    To-Do: implement smooth noise-robust differentiators for use on experimental data.
    http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
    """
    
    if x is None:
        x = np.arange(f.shape[axis]) * dx
    else:
        assert x.shape[0] == f.shape[axis]
    I = np.zeros(f.shape[axis])
    I[:2] = np.array([0, -1])
    I[-1] = 1
    I = circulant(I)
    I[0, 0] = -1
    I[-1, -1] = 1
    I[0, -1] = 0
    I[-1, 0] = 0
    H = np.zeros((f.shape[axis], 1))
    H[1:-1, 0] = x[2:] - x[:-2]
    H[0] = x[1] - x[0]
    H[-1] = x[-1] - x[-2]
    if axis == 0:
        return np.dot(I / H, f)
    else:
        return np.dot(I / H, f.T).T 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import circulant [as 别名]
def named_target_matrix(name, size):
    """
    Parameter:
        name: name of the target matrix
    Return:
        target_matrix: (n, n) numpy array for real matrices or (n, n, 2) for complex matrices.
    """
    if name == 'dft':
        return LA.dft(size, scale='sqrtn')[:, :, None].view('float64')
    elif name == 'idft':
        return np.ascontiguousarray(LA.dft(size, scale='sqrtn').conj().T)[:, :, None].view('float64')
    elif name == 'dft2':
        size_sr = int(math.sqrt(size))
        matrix = np.fft.fft2(np.eye(size_sr**2).reshape(-1, size_sr, size_sr), norm='ortho').reshape(-1, size_sr**2)
        # matrix1d = LA.dft(size_sr, scale='sqrtn')
        # assert np.allclose(np.kron(m1d, m1d), matrix)
        # return matrix[:, :, None].view('float64')
        from butterfly.utils import bitreversal_permutation
        br_perm = bitreversal_permutation(size_sr)
        br_perm2 = np.arange(size_sr**2).reshape(size_sr, size_sr)[br_perm][:, br_perm].reshape(-1)
        matrix = np.ascontiguousarray(matrix[:, br_perm2])
        return matrix[:, :, None].view('float64')
    elif name == 'dct':
        # Need to transpose as dct acts on rows of matrix np.eye, not columns
        # return dct(np.eye(size), norm='ortho').T
        return dct(np.eye(size)).T / math.sqrt(size)
    elif name == 'dst':
        return dst(np.eye(size)).T / math.sqrt(size)
    elif name == 'hadamard':
        return LA.hadamard(size) / math.sqrt(size)
    elif name == 'hadamard2':
        size_sr = int(math.sqrt(size))
        matrix1d = LA.hadamard(size_sr) / math.sqrt(size_sr)
        return np.kron(matrix1d, matrix1d)
    elif name == 'b2':
        size_sr = int(math.sqrt(size))
        import torch
        from butterfly import Block2x2DiagProduct
        b = Block2x2DiagProduct(size_sr)
        matrix1d = b(torch.eye(size_sr)).t().detach().numpy()
        return np.kron(matrix1d, matrix1d)
    elif name == 'convolution':
        np.random.seed(0)
        x = np.random.randn(size)
        return LA.circulant(x) / math.sqrt(size)
    elif name == 'hartley':
        return hartley_matrix(size) / math.sqrt(size)
    elif name == 'haar':
        return haar_matrix(size, normalized=True) / math.sqrt(size)
    elif name == 'legendre':
        grid = np.linspace(-1, 1, size + 2)[1:-1]
        return legendre.legvander(grid, size - 1).T / math.sqrt(size)
    elif name == 'hilbert':
        H = hilbert_matrix(size)
        return H / np.linalg.norm(H, 2)
    elif name == 'randn':
        np.random.seed(0)
        return np.random.randn(size, size) / math.sqrt(size)
    else:
        assert False, 'Target matrix name not recognized or implemented'