def test_radon_iradon():
    size = 100
    debug = False
    image = np.tri(size) + np.tri(size)[::-1]
    for filter_type in ["ramp", "shepp-logan", "cosine", "hamming", "hann"]:
        reconstructed = iradon(radon(image), filter=filter_type)

        image = rescale(image)
        reconstructed = rescale(reconstructed)
        delta = np.mean(np.abs(image - reconstructed))

        if debug:
            print(delta)
            import matplotlib.pyplot as plt
            f, (ax1, ax2) = plt.subplots(1, 2)
            ax1.imshow(image, cmap=plt.cm.gray)
            ax2.imshow(reconstructed, cmap=plt.cm.gray)
            plt.show()

        assert delta < 0.05

    reconstructed = iradon(radon(image), filter="ramp", interpolation="nearest")
    delta = np.mean(abs(image - reconstructed))
    assert delta < 0.05
    size = 20
    image = np.tri(size) + np.tri(size)[::-1]
    reconstructed = iradon(radon(image), filter="ramp", interpolation="nearest")

def diffusion_sm_shift_K(N,K,a_bar,a_degger_bar):
	# matrices to shift above and below diagonal and solve initial equation
	diag_n_count = diag_count(N)
	delta = np.eye(N)*a_bar - np.dot(diag_n_count,np.tri(N,N,-1)-np.tri(N,N,-2))
	delta_degger = np.eye(N)*a_degger_bar - np.tri(N,N,1)+np.tri(N,N,0)

	# store generating function coefficients in this matrix:
	gen = np.zeros((N,K))

	# two initial vectors for k2=0 and k2=1
	#for n1 in np.arange(N):
		#gen[n1,0] = 1./2.**(N-1)*np.float(math.factorial(N-1))/(math.factorial(N-1-n1)*math.factorial(n1))
	gen[:,0] = np.ones(N)
	gen[:,1] = np.dot(delta,gen[:,0])
	allcond = []
	# calculate full matrix gen, i.e. for row n1 and column k2 it contains gen(n1,k2) at steady state
	for k2 in np.arange(1,K-1):
		gen[:,k2+1] = la.solve((k2+1)*delta_degger , -np.dot( np.dot(delta_degger,delta) + np.eye(N)*k2 ,gen[:,k2]) - np.dot(delta,gen[:,k2-1]))
		eigval = la.eig((k2+1)*delta_degger)
		cond = np.float(np.max(np.abs(eigval[0])))/np.min(np.abs(eigval[0]))
#		print cond
		allcond.append(cond)
		
	coeff_mat = mat_K(a_bar*a_degger_bar,N,K).T

	p_mat = np.dot(gen,coeff_mat)
	p_n1 = np.sum(p_mat, axis=1) # sum of each row n gives value of the marginal distribution p(n)
	p_n2 = np.sum(p_mat, axis=0) # sum of each column m gives value of the marginal distribution p(m)

	n1_mean = np.dot(np.arange(N),p_n1) # compare with theoretically expected value: 
	n2_mean = np.dot(np.arange(N),p_n2)
	#print((n1_mean,n2_mean,np.sum(p_mat)))
	return (p_mat,allcond)

def test_iradon_angles():
    """
    Test with different number of projections
    """
    size = 100
    # Synthetic data
    image = np.tri(size) + np.tri(size)[::-1]
    # Large number of projections: a good quality is expected
    nb_angles = 200
    radon_image_200 = radon(image, theta=np.linspace(0, 180, nb_angles,
                                                     endpoint=False))
    reconstructed = iradon(radon_image_200)
    delta_200 = np.mean(abs(_rescale_intensity(image) - _rescale_intensity(reconstructed)))
    assert delta_200 < 0.03
    # Lower number of projections
    nb_angles = 80
    radon_image_80 = radon(image, theta=np.linspace(0, 180, nb_angles,
                                                    endpoint=False))
    # Test whether the sum of all projections is approximately the same
    s = radon_image_80.sum(axis=0)
    assert np.allclose(s, s[0], rtol=0.01)
    reconstructed = iradon(radon_image_80)
    delta_80 = np.mean(abs(image / np.max(image) -
                           reconstructed / np.max(reconstructed)))
    # Loss of quality when the number of projections is reduced
    assert delta_80 > delta_200

def tri_flat(array, UPPER=True):
    """
    Flattens the upper/lower triangle of a square matrix.
    
    Parameters
    ----------
    array : np.ndarray
        square matrix
        
    UPPER : boolean
        Upper or lower triangle to flatten (defaults to upper). If
        the matrix is symmetric, this parameter is unnecessary.
    
    Returns
    -------
    array : np.ndarray
        vector representation of the upper / lower triangle
    """
    
    C = array.shape[0]
    if UPPER:
        mask = np.asarray(np.invert(np.tri(C,C,dtype=bool)),dtype=float)
    else:
        mask = np.asarray(np.invert(np.tri(C,C,dtype=bool).transpose()),dtype=float)
        
    x,y = mask.nonzero()
    return array[x,y]

def diffusion_recurrence(N):
	diag = diag_count(N)
	A1 = np.tri(N,N,1)-np.tri(N,N,0)
	A2 = np.tri(N,N,2)-np.tri(N,N,1)

	p_mat_old = np.zeros((N,N))

	# two initial vectors:
	p_mat_old[:,0] = np.zeros(N)
	p_mat_old[N-1,0] = 1.
	p_mat_old[:,1] = np.dot(np.dot(A1,diag),p_mat_old[:,0])

	# calculate matrix containing the full distribution, i.e. for row n1 and column n2 it contains p(n1,n2) at steady state
	for n2 in np.arange(2,N):
		p_mat_old[:,n2] = 1./n2*(-np.dot(np.dot(A2,diag),p_mat_old[:,n2-2]) + (n2-1)*np.dot(A1,p_mat_old[:,n2-1]) + np.dot(np.dot(A1,diag),p_mat_old[:,n2-1]))

	p_mat_old = p_mat_old/np.sum(p_mat_old)
	return p_mat_old

#p_mat_old = diffusion_recurrence(N)

#p_n1_old = np.sum(p_mat_old, axis=1) # sum of each row n gives value of the marginal distribution p(n)
#p_n2_old = np.sum(p_mat_old, axis=0) # sum of each column m gives value of the marginal distribution p(m)

#n1_mean_old = np.dot(np.arange(N),p_n1_old) # compare with theoretically expected value: 
#n2_mean_old = np.dot(np.arange(N),p_n2_old)

def test_template():
    size = 100
    # Float prefactors ensure that image range is between 0 and 1
    image = np.full((400, 400), 0.5)
    target = 0.1 * (np.tri(size) + np.tri(size)[::-1])
    target_positions = [(50, 50), (200, 200)]
    for x, y in target_positions:
        image[x:x + size, y:y + size] = target
    np.random.seed(1)
    image += 0.1 * np.random.uniform(size=(400, 400))

    result = match_template(image, target)
    delta = 5

    positions = peak_local_max(result, min_distance=delta)

    if len(positions) > 2:
        # Keep the two maximum peaks.
        intensities = result[tuple(positions.T)]
        i_maxsort = np.argsort(intensities)[::-1]
        positions = positions[i_maxsort][:2]

    # Sort so that order matches `target_positions`.
    positions = positions[np.argsort(positions[:, 0])]

    for xy_target, xy in zip(target_positions, positions):
        assert_almost_equal(xy, xy_target)

def crankNicolson(condInitialesPhi, condInitialesPsi, condSpatiales = None, tMax = 0.001, dt=10**-6, v = 1, dx = 1):

	if np.size(condInitialesPhi) != np.size(condInitialesPsi) :
		raise Exception("La taille de condInitialesPhi doit être semblable à condInitialesPsi")

	# Constantes utiles
	n = np.size(condInitialesPhi)
	k = -dt * v**2 / dx**2 / 2
	N = int(tMax / dt)

	# Matrice de l’évolution du système 
	evolution = np.zeros((N+1,2*n))
	evolution[0,:n] = condInitialesPhi
	evolution[0,n:] = condInitialesPsi

	# On créer la matrice d'évolution 
	I = np.eye(n)
	A = np.tri(n, k = 1).T * np.tri(n, k=-1)
	A = (A + A.T - 2 * I) * k 
	M = np.array(np.bmat(((I, -dt*I/2),(A, I))))
	K = np.array(np.bmat(((I, dt*I/2),(-A, I))))
	invM = np.linalg.inv(M)
	matriceEvolution = np.dot(invM,K)	

	# On applique les conditions spatiales obtenant la liste des points qui varie dans le temps.
	if condSpatiales is not None :
		matriceEvolution[condSpatiales] = np.zeros(2*n)
		matriceEvolution[condSpatiales, condSpatiales] = 1
		matriceEvolution[condSpatiales+n] = np.zeros(2*n)

	for i in range(1,N+1):
		evolution[i] = np.dot(matriceEvolution,evolution[i-1])

	return evolution[:,:n], evolution[:,n:]

def get_potentials_wake(ring, harm_cav, lamb):
    E0 = ring.E0
    C0 = ring.C0
    T0 = ring.T0
    h = ring.harm_num
    Rs = harm_cav.shunt_imp
    alpha = harm_cav.alpha
    beta = harm_cav.beta
    z = lamb.z
    dist = lamb.dist*lamb.cur[:, None]

    Ll = 1/(1-np.exp(-beta*C0))  # Lesser
    Gl = Ll*np.exp(-beta*C0)     # Greater or Equal

    A = Ll*np.tri(h, h, -1) + Gl*np.tri(h, h).T
    ind_b = np.arange(h)
    dif = ind_b[:, None] - ind_b[None, :]
    B = np.exp(-beta*C0*dif/h)
    lamb_w = np.trapz(dist*np.exp(beta*z)[None, :], z)
    V = np.dot(A*B, lamb_w)

    Sn = scy_int.cumtrapz(np.exp(beta*z)[None, :] * dist, x=z, initial=0*1j)
    Vt = np.exp(-beta*z)[None, :] * (V[:, None] + Sn)
    k = 1

    return -T0/E0 * 2*alpha*Rs*(Vt.real - alpha/harm_cav.wrb*Vt.imag), k

def test_roberts_diagonal1():
    """Roberts' on an edge should be a one diagonal"""
    image = np.tri(10, 10, 0)
    expected = ~(np.tri(10, 10, -1).astype(bool) | \
                 np.tri(10, 10, -2).astype(bool).transpose())
    expected = _mask_filter_result(expected,None)
    result = F.roberts(image).astype(bool)
    assert_close(result,expected)

def test_roberts_diagonal2():
    """Roberts' filter on a diagonal edge should be a diagonal line."""
    image = np.rot90(np.tri(10, 10, 0), 3)
    expected = ~np.rot90(np.tri(10, 10, -1).astype(bool) |
                         np.tri(10, 10, -2).astype(bool).transpose())
    expected = _mask_filter_result(expected, None)
    result = F.roberts(image).astype(bool)
    assert_close(result, expected)

def test_roberts_diagonal1():
    """Roberts' filter on a diagonal edge should be a diagonal line."""
    image = np.tri(10, 10, 0)
    expected = ~(np.tri(10, 10, -1).astype(bool) |
                 np.tri(10, 10, -2).astype(bool).transpose())
    expected = _mask_filter_result(expected, None)
    result = filters.roberts(image).astype(bool)
    assert_array_almost_equal(result, expected)

def _make_cubic_spline_matrix(nn):
    ii = np.eye(nn,nn)
    sub_ii = np.tri(nn,nn,k=-1) - np.tri(nn,nn,k=-2)

    coeff = 4*ii + sub_ii + sub_ii.transpose()

    coeff[0][0] = 2
    coeff[nn-1][nn-1] = 2

    return np.matrix(coeff)

def test_roberts_diagonal2():
    """Roberts' on an edge should be a other diagonal"""
    diagonal = np.tri(10, 10, 0,dtype=int)
    rev_diagonal = np.rot90(diagonal.transpose(),1)

    image = (rev_diagonal > 0).astype(float)
    expected = ~np.rot90((np.tri(10, 10, -1).astype(bool) | \
                np.tri(10, 10, -2).astype(bool).transpose()),1)
    expected = _mask_filter_result(expected,None)
    result = F.roberts(image).astype(bool)
    assert_close(result,expected)

def subsequence_ecdf(x, w):
    n = x.shape[0]
    return numpy.sum(
        numpy.where(
            (1 - numpy.tri(n, k=-1, dtype=bool)).reshape((n, 1, n)).repeat(n, axis=1)
            & numpy.tri(n, dtype=bool).reshape((1, n, n)).repeat(n, axis=0)
            & (x.reshape((n, 1)).repeat(n, axis=1) >= x).reshape((1, n, n)).repeat(n, axis=0),
            w,
            0,
        ),
        axis=2,
    ) / subsequence_sums(w)

def triDiag(n):
	"""Retourne une matrice n x n de la forme 
		A = [-2 1 0 0 ...]
			[1 -2 1 0 ...]
			[0 1 -2 1 ...]
			[ .......... ]
			[... 0 1 -2 1]
			[... 0 0 1 -2]
	"""
	A = np.tri(n, k = 1).T * np.tri(n, k=-1)
	A = (A + A.T - 2 * np.eye(n))
	return A

def test_frt():
    SIZE = 59
    try:
        import sympy.ntheory as sn
        assert sn.isprime(SIZE) == True
    except ImportError:
        pass

    # Generate a test image
    L = np.tri(SIZE, dtype=np.int32) + np.tri(SIZE, dtype=np.int32)[::-1]
    f = frt2(L)
    fi = ifrt2(f)
    assert len(np.nonzero(L - fi)[0]) == 0

def get_gauss_means_vec(u1, v1, I):
    n = len(I)
    aux = np.repeat(u1[:,np.newaxis],n, axis=1)

    #calculate x*\exp{-\int_s^t g(v)dv}
    aux[np.tri(n, n, 0)==0] = 1 
    mu1 = np.cumprod(aux, 0)

    aux[np.tri(n, n, -1)==0] = 1 
    y = np.cumprod(aux,0)
    y[np.tri(n, n, 0)==0] = 0
    y = v1*I*y
    mu2 = np.cumsum(y[:,::-1], 1)[:, ::-1]
    return mu1, mu2

def vech(X, order='F'):
    """Returns vectorization of lower triangle of X in column major order by default."""
    assert X.ndim == 2, 'vech operator requires a matrix.'
    
    m, n = X.shape
    if order == 'F':
        idx = np.where(1 - np.tri(n,m, -1, dtype=int))
        return X.T[idx]
    elif order == 'C':
        i,j = np.where(np.tri(m,n, dtype=int))
    else:
        raise Exception("Only order C and F are allowed")
    
    return X[i,j]

def crankNicolson2D(condInitialesPhi, condInitialesPsi, condSpatiales = None, tMax = 0.001, dt=10**-6, v = 100, dx = 1, intervalSauvegarde=1):
	if not np.array_equal(np.shape(condInitialesPhi), np.shape(condInitialesPsi)):
		raise Exception("La taille de condInitialesPhi doit être similaire à la taille condInitialesPsi.")


	# Constantes utiles
	n = np.shape(condInitialesPhi)  #(ny,nx)

	if n[0] != n[1] :
		raise Exception("Les dimensions x et y doivent être similaires.")

	n = n[0]
	k = -dt * v**2 / dx**2 / 2
	N = int(tMax / dt)

	# Matrice de l’évolution du système 
	evolution = np.zeros((int(N/intervalSauvegarde)+1,2*n*n))
	evolution[0,:n*n] = condInitialesPhi.flatten()
	evolution[0,n*n:] = condInitialesPsi.flatten()

	phi = evolution[0]

	I = np.eye(n*n)
	A = np.tri(n*n, k = 1).T * np.tri(n*n, k=-1)
	A = (A + A.T - 4 * np.eye(n*n))*k
	B = np.eye((n-1)*n)*k
	A[:-n,n:] = A[:-n,n:] + B
	A[n:,:-n] = A[n:,:-n] + B

	M = np.array(np.bmat( ((I, -dt*I/2), (A, I))))
	K = np.array(np.bmat( ((I, dt*I/2), (-A, I))))

	invM = np.linalg.inv(M)
	matriceEvolution = np.dot(invM,K)	

	# On applique les conditions spatiales obtenant la liste des points qui varie dans le temps.
	if condSpatiales is not None :
		idx = np.array(condSpatiales[0]+n*condSpatiales[1],"int")
		matriceEvolution[idx] = np.zeros(2*n*n)
		matriceEvolution[idx, idx] = 1
		matriceEvolution[idx+n*n] = np.zeros(2*n*n)

	for i in range(1,N+1):
		phi = np.dot(matriceEvolution,phi)
		if i % intervalSauvegarde == 0:
			evolution[int(i // intervalSauvegarde)] = phi

	return evolution[:,:n*n], evolution[:,n*n:]

def triang_decomp(c):
    """Return a lower triangular matrix B that B * B.T = C"""
    n = c.shape[0]
    b = np.tri(n, n)

    b[1:n, 0] = c[1:n, 0]

    for i in xrange(1, n):
        b[i, 1:i+1] = np.sqrt(1 - c[i, 0]**2)

    for i in xrange(2, n):
        for j in xrange(1, i):
            b1 = np.dot(b[j, 0:j], b[i, 0:j].T)
            b2 = np.dot(b[j, j], b[i, j])
            cosinv = (c[i, j] - b1)/b2

            if np.isfinite(cosinv):
                if cosinv > 1:
                    b[i, j] = b[i, j]
                    b[i, j+1:n+1] = 0
                elif cosinv < -1:
                    b[i, j] = -b[i, j]
                    b[i, j+1:n+1] = 0
                else:
                    b[i, j] = b[i, j]*cosinv
                    sinTheta = np.sqrt(1 - cosinv**2)
                    for k in xrange(j+1, n):
                        b[i, k] = b[i, k]*sinTheta

    return b