Python scipy.linalg.circulant() Examples

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 

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 

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 

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 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 

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)) 

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 

def _setup(self, config):
        torch.manual_seed(config['seed'])
        self.model = ButterflyProduct(size=config['size'],
                                      complex=True,
                                      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']
        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' 

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' 

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)) 

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)) 

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) ) 

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)) 

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) 

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' 

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 

def jl_fit_transform(X, objective_dim, method="basic"):
    """Fit and transform the input data by Johnson–Lindenstrauss process.
    See :cite:`johnson1984extensions` for details.

    Parameters
    ----------
    X : numpy array of shape (n_samples, n_features)
        The input samples.

    objective_dim : int
        The expected output dimension.

    method : string, optional (default = 'basic')
        The JL projection method:

        - "basic": each component of the transformation matrix is taken at
          random in N(0,1).
        - "discrete", each component of the transformation matrix is taken at
          random in {-1,1}.
        - "circulant": the first row of the transformation matrix is taken at
          random in N(0,1), and each row is obtained from the previous one
          by a one-left shift.
        - "toeplitz": the first row and column of the transformation matrix
          is taken at random in N(0,1), and each diagonal has a constant value
          taken from these first vector.

    Returns
    -------
    X_transformed : numpy array of shape (n_samples, objective_dim)
        The dataset after the JL projection.

    jl_transformer : object
        Transformer instance.
    """
    if method.lower() == "basic":
        jl_transformer = (1 / math.sqrt(objective_dim)) \
                          * np.random.normal(0, 1, size=(objective_dim,
                                                         len(X[0])))
    elif method.lower() == "discrete":
        jl_transformer = (1 / math.sqrt(objective_dim)) \
                          * np.random.choice([-1, 1],
                                             size=(objective_dim, len(X[0])))
    elif method.lower() == "circulant":
        from scipy.linalg import circulant
        first_row = np.random.normal(0, 1, size=(1, len(X[0])))
        jl_transformer = ((1 / math.sqrt(objective_dim))
                           * circulant(first_row))[:objective_dim]
    elif method.lower() == "toeplitz":
        from scipy.linalg import toeplitz
        first_row = np.random.normal(0, 1, size=(1, len(X[0])))
        first_column = np.random.normal(0, 1, size=(1, objective_dim))
        jl_transformer = ((1 / math.sqrt(objective_dim))
                           * toeplitz(first_column, first_row))
    else:
        NotImplementedError('Wrong transformation type')

    jl_transformer = jl_transformer.T

    return np.dot(X, jl_transformer), jl_transformer 

def gram_block_circulant(ds, n_times_valid, method='full',
                         sample_weights=None):
    """Returns ...

    Parameters
    ----------
    ds : array, shape (n_atoms, n_times_atom)
        The atoms
    n_times_valid : int
        n_times - n_times_atom + 1
    method : string
        If 'full', returns full circulant matrix.
        If 'scipy', returns scipy linear operator.
        If 'custom', returns custom linear operator.
    sample_weights : array, shape (n_times, )
        The sample weights for one trial
    """
    from scipy.sparse.linalg import LinearOperator
    from functools import partial

    n_atoms, n_times_atom = ds.shape
    n_times = n_times_valid + n_times_atom - 1

    if method == 'full':
        D = np.zeros((n_times, n_atoms * n_times_valid))
        for k_idx in range(n_atoms):
            d_padded = np.zeros((n_times, ))
            d_padded[:n_times_atom] = ds[k_idx]
            start = k_idx * n_times_valid
            stop = start + n_times_valid
            D[:, start:stop] = linalg.circulant((d_padded))[:, :n_times_valid]
        if sample_weights is not None:
            wD = sample_weights[:, None] * D
            return np.dot(D.T, wD)
        else:
            return np.dot(D.T, D)

    elif method == 'scipy':

        def matvec(v, ds):
            assert v.shape[0] % ds.shape[0] == 0
            return _fprime(ds, v, Xi=None, reg=None,
                           sample_weights=sample_weights)

        D = LinearOperator((n_atoms * n_times_valid, n_atoms * n_times_valid),
                           matvec=partial(matvec, ds=ds))

    elif method == 'custom':
        return CustomLinearOperator(ds, n_times_valid, sample_weights)
    else:
        raise ValueError('Unkown method %s.' % method)
    return D 

def block_circ(a):
    '''generate full representation of 2-level circulant matrix

    Args:
        a (ndarray): 1-st column of the circulant matrix in proper shape.
    Returns:
        Full filled circulant 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)
            
    if a.ndim == 1:
        return splin.circulant(a)
    n_total = np.prod(a.shape)
    a_shape = a.shape
    A = np.zeros(np.hstack([np.array(a.shape),np.array(a.shape)]))
    A_shape = A.shape
    x_slice = [0]*a.ndim
    x_slice[-1] = slice(None)
    x_target_slice = [0]*a.ndim
    x_target_slice[-1] = slice(None)
    y_slice = [slice(None)]*2
    a = a.reshape(-1,a_shape[-1])
    A = A.reshape(-1,a_shape[-1],*a_shape)
    for i,sub_column in enumerate(a):
        print(sub_column)
        y_slice[0] = i
        A[tuple(y_slice+x_slice)] = splin.circulant(sub_column)
    for ilevel in range(len(a_shape)-1):
        A = A.reshape(-1,*a_shape[len(a_shape)-ilevel-2:],*a_shape)
        y_slice = [slice(None)]*(ilevel+3)
        y_target_slice = [slice(None)]*(ilevel+3)
        for k in range(A.shape[0]):
            y_slice[0] = k
            y_target_slice[0] = k
            for i in range(a_shape[len(a_shape)-ilevel-2]):
                y_slice[1] = i
                for j in range(1,a_shape[len(a_shape)-ilevel-2]):
                    x_slice[len(a_shape)-ilevel-2] = j
                    y_target_slice[1] = np.mod(i-j,a_shape[len(a_shape)-ilevel-2])
                    A[tuple(y_slice+x_slice)] = A[tuple(y_target_slice+x_target_slice)]
        x_slice[len(a_shape)-ilevel-2] = slice(None)
        x_target_slice[len(a_shape)-ilevel-2] =slice(None)
    A = A.reshape(A_shape)
    return A