Python linalg.circulant函数代码示例

def display(n=10):
    s1 = [n-i for i in range(n)]+[0 for i in range(n-1)]
    # build causal right circulant matrix
    # np.matrix rather than np.ndarray could perform matrix mulitplication
    s1 = np.matrix(np.flipud(np.fliplr(circulant(s1)[n-1:,:])))
    #s2 = np.matrix(s1).T
    # using a diffrent and random matrix for generality
    s2 = [np.random.randint(1,20) for i in range(n)]+[0 for i in range(n-1)]
    s2 = np.matrix(np.flipud(np.fliplr(circulant(s2)[n-1:,:]))).T

    print s1
    print "\n", s2
    print "\n", s1*s2

def spectral_diff_matrix(N,dt,diff):
  ''' 
  generates a periodic sinc differentation matrix. This is equivalent 
  
  Parameters
  ----------
    N : number of observations 
    dt : sample spacing
    diff : derivative order (max=2)

  '''
  scale = dt*N/(2*np.pi)
  dt = 2*np.pi/N
  t,h = sympy.symbols('t,h')
  sinc = sympy.sin(sympy.pi*t/h)/((2*sympy.pi/h)*sympy.tan(t/2))
  if diff == 0:
    sinc_diff = sinc
  else:
    sinc_diff = sinc.diff(*(t,)*diff)

  func = sympy.lambdify((t,h),sinc_diff,'numpy')
  times = dt*np.arange(N)
  val = func(times,dt)
  if diff == 0:
    val[0] = 1.0
  elif diff == 1:
    val[0] = 0.0
  elif diff == 2:
    val[0] = -(np.pi**2/(3*dt**2)) - 1.0/6.0

  D = circulant(val)/scale**diff
  return D

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 _distance_matrix(L):
    Dmax = L//2
 
    D  = range(Dmax+1)
    D += D[-2+(L%2):0:-1]
 
    return circulant(D)/Dmax

 def test_basic3(self):
     # b is a 3-d matrix.
     c = np.array([1, 2, -3, -5])
     b = np.arange(24).reshape(4, 3, 2)
     x = solve_circulant(c, b)
     y = solve(circulant(c), b)
     assert_allclose(x, y)

 def test_complex(self):
     # Complex b and c
     c = np.array([1+2j, -3, 4j, 5])
     b = np.arange(8).reshape(4, 2) + 0.5j
     x = solve_circulant(c, b)
     y = solve(circulant(c), b)
     assert_allclose(x, y)

def test_toeplitz(N=50):
    print("Testing circulant linear algebra...")
    x = np.linspace(0, 10, N)
    y = np.vstack((np.sin(x), np.cos(x), x, x**2)).T
    c_row = np.exp(-0.5 * x ** 2)
    c_row[0] += 0.1
    cnum = circulant(c_row)
    cmat = CirculantMatrix(c_row)

    # Test dot products.
    assert np.allclose(np.dot(cnum, y[:, 0]), cmat.dot(y[:, 0]))
    assert np.allclose(np.dot(cnum, y), cmat.dot(y))

    # Test solves.
    assert np.allclose(np.linalg.solve(cnum, y[:, 0]), cmat.solve(y[:, 0]))
    assert np.allclose(np.linalg.solve(cnum, y), cmat.solve(y))

    # Test eigenvalues.
    ev = np.linalg.eigvals(cnum)
    ev = ev[np.argsort(np.abs(ev))[::-1]]
    assert np.allclose(np.abs(cmat.eigvals()), np.abs(ev))

    print("Testing Toeplitz linear algebra...")
    tnum = toeplitz(c_row)
    tmat = ToeplitzMatrix(c_row)

    # Test dot products.
    assert np.allclose(np.dot(tnum, y[:, 0]), tmat.dot(y[:, 0]))
    assert np.allclose(np.dot(tnum, y), tmat.dot(y))

    # Test solves.
    assert np.allclose(np.linalg.solve(tnum, y[:, 0]),
                       tmat.solve(y[:, 0], tol=1e-12, verbose=True))
    assert np.allclose(np.linalg.solve(tnum, y),
                       tmat.solve(y, tol=1e-12, verbose=True))

def calculation(data, m, n, activation, w):
    nrow, ncol = sorted((m, n))
    r = np.random.rand(nrow, 1)
    circul_matrix = circulant(r)
    circul_matrix = np.hstack((circul_matrix, circul_matrix[:, :ncol-nrow]))

    fft_r = np.fft.fft(r)
    fft_x = np.fft.fft(data)
    Rx = np.fft.ifft(fft_r * fft_x)
    if activation == 1:
        hx = sigmoid(Rx)
    elif activation == 2:
        hx = np.tanh(Rx)

    FP = hx
    rev_x = np.flipud(data)
    s_rev_x = np.roll(rev_x, 1, axis=0)
    if activation == 1:
        dhx = hx * (1 - hx)
    elif activation == 2:
        dhx = 1 - hx**2

    fft_s_rev_x = np.fft.fft(s_rev_x.transpose())
    fft_wT_rox = np.fft.fft((w * dhx).transpose())
    BP = np.fft.ifft((fft_s_rev_x * fft_wT_rox).transpose())

    return circul_matrix, FP, BP

def est_covs(samples, epsilon):
    c1 = samples[:,:1000]
    cf = samples[:,-1000:]
    
    r1 = np.array([2+epsilon, -1, 0, -1])
    Q = sp.circulant(r1)

    return 4-np.trace(np.dot(Q,np.cov(c1))), 4-np.trace(np.dot(Q, np.cov(cf))), 4-np.trace(np.dot(Q,np.cov(samples)))

def generate_all_shifts(atom, signal_atom_diff):
    """ Shifted version of a matrix with the number
    of shifts being signal_atom_diff
    the atom is a vector and the shifted versions 
    are the rows
    extra zeros are tacked on
    """
    return circulant(np.hstack((atom, np.zeros(signal_atom_diff)))).T[:signal_atom_diff+1]

 def test_random_b_and_c(self):
     # Random b and c
     np.random.seed(54321)
     c = np.random.randn(50)
     b = np.random.randn(50)
     x = solve_circulant(c, b)
     y = solve(circulant(c), b)
     assert_allclose(x, y)

def C(size):
    
    firstrow = np.zeros(size)
    firstrow[-1] = -1
    firstrow[0] = 2
    firstrow[1] = -1
    
    return linalg.circulant(firstrow)

def Dc(size): #periodic first derivative
    
    firstrow = np.zeros(size)
    firstrow[-1] = -1
    firstrow[0] = 1

    
    return linalg.circulant(firstrow)

 def test_singular(self):
     # c gives a singular circulant matrix.
     c = np.array([1, 1, 0, 0])
     b = np.array([1, 2, 3, 4])
     x = solve_circulant(c, b, singular='lstsq')
     y, res, rnk, s = lstsq(circulant(c), b)
     assert_allclose(x, y)
     assert_raises(LinAlgError, solve_circulant, x, y)

def make_gibbs_hist(epsilon):
    total_rhos = np.zeros(100)

    r1 = np.array([2+epsilon, -1, 0, -1])
    Q = sp.circulant(r1)
    for i in range(100):
        sample = run_gibbs_normal(10000, epsilon)
        total_rhos[i] = 4-np.trace(np.dot(Q, np.cov(sample)))
    plt.hist(total_rhos)