Python numpy.tri函数代码示例

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