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)

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

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

 def prepare_autoregression_coefficients(self,CF):
     '''Calculate autoregression coefficients (cVec) and variance of 
     additional noise term (sigma**2).
     Is it possible to simplify the matrix inverse in case of a circulant
     matrix?'''
     B = circulant(CF[:-1]) # up to k-1
     Binv = np.linalg.inv(B) # invert matrix
     cVec = np.dot(Binv,CF[1:]) # from 1 to k
     sigma = np.sqrt(CF[0] - np.dot(cVec,CF[1:]))
     
     return sigma, cVec

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 __init__(self, a, nu, alpha, v0, xaxis):
   self.a     = a
   self.nu    = nu
   self.alpha = alpha
   self.v0    = v0
   self.xaxis = xaxis
   self.dim   = np.size(xaxis)
   self.dx    = xaxis[1] - xaxis[0]
   e     = np.zeros(self.dim)
   e[0]  = -2.0
   e[1]  = 1.0
   e[-1] = 1.0
   self.S  = sp.csc_matrix(spla.circulant(e))
   self.S *= (self.nu/self.dx**2)

    def test_axis_args(self):
        # Test use of caxis, baxis and outaxis.

        # c has shape (2, 1, 4)
        c = np.array([[[-1, 2.5, 3, 3.5]], [[1, 6, 6, 6.5]]])

        # b has shape (3, 4)
        b = np.array([[0, 0, 1, 1], [1, 1, 0, 0], [1, -1, 0, 0]])

        x = solve_circulant(c, b, baxis=1)
        assert_equal(x.shape, (4, 2, 3))
        expected = np.empty_like(x)
        expected[:, 0, :] = solve(circulant(c[0]), b.T)
        expected[:, 1, :] = solve(circulant(c[1]), b.T)
        assert_allclose(x, expected)

        x = solve_circulant(c, b, baxis=1, outaxis=-1)
        assert_equal(x.shape, (2, 3, 4))
        assert_allclose(np.rollaxis(x, -1), expected)

        # np.swapaxes(c, 1, 2) has shape (2, 4, 1); b.T has shape (4, 3).
        x = solve_circulant(np.swapaxes(c, 1, 2), b.T, caxis=1)
        assert_equal(x.shape, (4, 2, 3))
        assert_allclose(x, expected)