Python numpy.triu函数代码示例

def vRaman(x,omega=1.0,delta=0.0,epsilon=0.048,phi=4.0/3.0):
    x=np.array(x)    
    s=1
    v=v=np.einsum('i,jk->ijk',omega*np.exp(1.0j*2.0*phi*x),Fx(s)*np.sqrt(2.0)/2.0)
    v=np.array([np.triu(v[i])+np.conjugate(np.triu(v[i],1)).transpose() for i in range(x.size)])
    v+=np.array([np.diag([epsilon+delta,0.0,epsilon-delta])]*x.size)
    return v

 def pressure(self):
     ps = 1.
     sig = 1.
     eps = 1.
     d = 2.
     N = np.size(self.c.x)
     
     dx = self.c.dx()
     dy = self.c.dy()
     dz = self.c.dz()
     dr2 = self.c.dr2()
     
     r_mag = (dx ** 2 + dy ** 2 + dz ** 2)
     r_mag = np.nan_to_num(r_mag)
     
     px = dx * (24 * eps / r_mag) * ((2 * (sig**12 / r_mag**6)) - (sig**6 / r_mag**3))
     px = np.nan_to_num(px)
     px = np.triu(px)
     py = dy * (24 * eps / r_mag) * ((2 * (sig**12 / r_mag**6)) - (sig**6 / r_mag**3))
     py = np.nan_to_num(py)
     py = np.triu(py)
     pz = dz * (24 * eps / r_mag) * ((2 * (sig**12 / r_mag**6)) - (sig**6 / r_mag**3))
     pz = np.nan_to_num(pz)
     pz = np.triu(pz)
     pt = px + py + pz
     return 1 / (d * N * self.ke()) * np.sum(pt)

    def test_dynamic_programming_logic(self):
        # Test for the dynamic programming part
        # This test is directly taken from Cormen page 376.
        arrays = [np.random.random((30, 35)),
                  np.random.random((35, 15)),
                  np.random.random((15, 5)),
                  np.random.random((5, 10)),
                  np.random.random((10, 20)),
                  np.random.random((20, 25))]
        m_expected = np.array([[0., 15750., 7875., 9375., 11875., 15125.],
                               [0.,     0., 2625., 4375.,  7125., 10500.],
                               [0.,     0.,    0.,  750.,  2500.,  5375.],
                               [0.,     0.,    0.,    0.,  1000.,  3500.],
                               [0.,     0.,    0.,    0.,     0.,  5000.],
                               [0.,     0.,    0.,    0.,     0.,     0.]])
        s_expected = np.array([[0,  1,  1,  3,  3,  3],
                               [0,  0,  2,  3,  3,  3],
                               [0,  0,  0,  3,  3,  3],
                               [0,  0,  0,  0,  4,  5],
                               [0,  0,  0,  0,  0,  5],
                               [0,  0,  0,  0,  0,  0]], dtype=np.int)
        s_expected -= 1  # Cormen uses 1-based index, python does not.

        s, m = _multi_dot_matrix_chain_order(arrays, return_costs=True)

        # Only the upper triangular part (without the diagonal) is interesting.
        assert_almost_equal(np.triu(s[:-1, 1:]),
                            np.triu(s_expected[:-1, 1:]))
        assert_almost_equal(np.triu(m), np.triu(m_expected))

    def get_representational_similarity(self, corr='spearman'):
        logging.debug("Calculating representational similarity")

        response = self.response[:, 1:, :self.numbercells, 0] # orientation x freq x phase x cell, no blank
        response = response.reshape(self.number_ori * (self.number_tf-1), self.numbercells)
        Nstim, N = response.shape

        rep_sim = np.zeros((Nstim, Nstim))
        rep_sim_p = np.empty((Nstim, Nstim))
        if corr == 'pearson':
            for i in range(Nstim):
                for j in range(i, Nstim): # matrix is symmetric
                    rep_sim[i, j], rep_sim_p[i, j] = st.pearsonr(response[i], response[j])

        elif corr == 'spearman':
            for i in range(Nstim):
                for j in range(i, Nstim): # matrix is symmetric
                    rep_sim[i, j], rep_sim_p[i, j] = st.spearmanr(response[i], response[j])

        else:
            raise Exception('correlation should be pearson or spearman')

        rep_sim = np.triu(rep_sim) + np.triu(rep_sim, 1).T # fill in lower triangle
        rep_sim_p = np.triu(rep_sim_p) + np.triu(rep_sim_p, 1).T  # fill in lower triangle

        return rep_sim, rep_sim_p

    def test_tpsv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 10
            x = rand(n).astype(dtype)
            # Upper triangular array
            A = triu(rand(n, n)) if ind < 2 else triu(rand(n, n)+rand(n, n)*1j)
            A += eye(n)
            # Form the packed storage
            c, r = tril_indices(n)
            Ap = A[r, c]
            func, = get_blas_funcs(('tpsv',), dtype=dtype)

            y1 = func(n=n, ap=Ap, x=x)
            y2 = solve(A, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1)
            A[arange(n), arange(n)] = dtype(1)
            y2 = solve(A, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1, trans=1)
            y2 = solve(A.T, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1, trans=2)
            y2 = solve(A.conj().T, x)
            assert_array_almost_equal(y1, y2)

    def __init__(self, endmembers, alphas, energy_interaction, volume_interaction=None, entropy_interaction=None):

        self.n_endmembers = len(endmembers)

        # Create array of van Laar parameters
        self.alphas = np.array(alphas)

        # Create 2D arrays of interaction parameters
        self.We = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
        self.We[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in energy_interaction
                                                                for i in row])

        if entropy_interaction is not None:
            self.Ws = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
            self.Ws[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in entropy_interaction
                                                                        for i in row])
        else:
            self.Ws = np.zeros((self.n_endmembers, self.n_endmembers))

        if volume_interaction is not None:
            self.Wv = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
            self.Wv[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in volume_interaction
                                                                        for i in row])
        else:
            self.Wv = np.zeros((self.n_endmembers, self.n_endmembers))


        # initialize ideal solution model
        IdealSolution.__init__(self, endmembers)

    def compute_distances_and_pairs(self, pdb_file, nr_contacts=None, nr_noncontacts=None):
        #distance and contacts
        self.features['pair']['Cbdist'] = pdb.distance_map(pdb_file, self.L)

        #mask positions that have too many gaps
        gap_freq = 1 - (self.Ni / self.neff)
        highly_gapped_pos = np.where(gap_freq > self.max_gap_percentage)[0]
        self.features['pair']['Cbdist'][:,highly_gapped_pos] = np.nan
        self.features['pair']['Cbdist'][highly_gapped_pos, :] = np.nan

        #if there are unresolved residues, there will be nan in the distance_map
        with np.errstate(invalid='ignore'):
            self.features['pair']['contact'] = (self.features['pair']['Cbdist'] <= self.contact_threshold) * 1
            self.features['pair']['nocontact'] = (self.features['pair']['Cbdist'] > self.non_contact_threshold) * 1

        indices_contact = np.where(np.triu(self.features['pair']['contact'], k=self.seq_separation))
        indices_contact = tuple(shuffle(indices_contact[0],indices_contact[1], random_state=0))
        if nr_contacts:
            indices_contact = indices_contact[0][:nr_contacts], indices_contact[1][:nr_contacts]

        indices_nocontact = np.where(np.triu(self.features['pair']['nocontact'], k=self.seq_separation))
        indices_nocontact = tuple(shuffle(indices_nocontact[0],indices_nocontact[1], random_state=0))
        if nr_noncontacts:
            indices_nocontact = indices_nocontact[0][:nr_noncontacts], indices_nocontact[1][:nr_noncontacts]


        #update indices of i<j for only relevant pairs
        self.ij_ind_upper = np.array(list(indices_contact[0]) + list(indices_nocontact[0])), np.array(list(indices_contact[1]) + list(indices_nocontact[1]))

	def __init__(self, num_visibles, num_hiddens):
		"""
		Initializes the parameters of the SemiRBM.

		@type  num_visibles: integer
		@param num_visibles: number of visible units
		@type  num_hiddens:  integer
		@param num_hiddens:  number of hidden units
		"""

		AbstractBM.__init__(self, num_visibles, num_hiddens)

		# additional hyperparameters
		self.learning_rate_lateral = 0.01
		self.momentum_lateral = 0.5
		self.weight_decay_lateral = 0.

		self.damping = 0.2
		self.num_lateral_updates = 20
		self.sampling_method = AbstractBM.MF

		# additional parameters
		self.L = np.matrix(np.random.randn(num_visibles, num_visibles)) / num_visibles / 200
		self.L = np.triu(self.L) + np.triu(self.L).T - 2. * np.diag(np.diag(self.L))
		self.dL = np.zeros_like(self.L)

def calculate2_pseudoV(pred, truth, rnd=0.01, full_matrix=True, sym=False):
    if full_matrix:
        pred_cp = pred
        truth_cp = truth
    else: # make matrix upper triangular
        pred_cp = np.triu(pred)
        truth_cp = np.triu(truth)

    # Avoid dividing by zero by rounding everything less than rnd up to rnd
    # Note: it is ok to do this after making the matrix upper triangular
    # since the bottom triangle of the matrix will not affect the score

    size = np.array(pred_cp.shape)[1]
    res = 0 # result to be returned

    # do one row at a time to reduce memory usage
    for x in xrange(size):
        # (1 - rnd) will cast the pred_cp/truth_cp matrices automatically if they are int8
        pred_row = (1 - rnd) * pred_cp[x, ] + rnd
        truth_row = (1 - rnd) * truth_cp[x, ] + rnd

        pred_row /= np.sum(pred_row)
        truth_row /= np.sum(truth_row)
        if sym:
            res += np.sum(truth_row * np.log(truth_row/pred_row)) + np.sum(pred_row * np.log(pred_row/truth_row))
        else:
            res += np.sum(truth_row * np.log(truth_row/pred_row))
    return res

def vRaman(x,omega=1.0,delta=0.0,epsilon=0.048,phi=4.0/3.0):
    x=np.array(x)    
    s=1
    v=np.outer(omega*np.exp(1.0j*2.0*phi*x),Fx(s)*np.sqrt(2.0)/2.0).reshape(x.size,2*s+1,2*s+1)
    v=sLA.block_diag(*[np.triu(v[i])+np.conjugate(np.triu(v[i],1)).transpose() for i in range(x.size)])
    v+=sLA.block_diag(*[np.diag([epsilon-delta,0.0,epsilon+delta])]*x.size)
    return v

 def test_compute_modes(self):
     ws = N.identity(self.num_states)
     tol = 1e-6
     weights_full = N.mat(N.random.random((self.num_states, self.num_states)))
     weights_full = N.triu(weights_full) + N.triu(weights_full, 1).H
     weights_full = weights_full*weights_full
     weights_diag = N.random.random(self.num_states)
     weights_list = [None, weights_diag, weights_full]
     vec_array = N.random.random((self.num_states, self.num_vecs))
     for weights in weights_list:
         IP = VectorSpaceMatrices(weights=weights).compute_inner_product_mat
         correlation_mat_true = IP(vec_array, vec_array)
         eigen_vals_true, eigen_vecs_true = util.eigh(correlation_mat_true)
         build_coeff_mat_true = eigen_vecs_true * N.mat(N.diag(
             eigen_vals_true**-0.5))
         modes_true = vec_array.dot(build_coeff_mat_true)
         
         modes, eigen_vals, eigen_vecs, correlation_mat = \
             compute_POD_matrices_snaps_method(vec_array, self.mode_indices, 
             inner_product_weights=weights, return_all=True)
         
         N.testing.assert_allclose(eigen_vals, eigen_vals_true, rtol=tol)
         N.testing.assert_allclose(eigen_vecs, eigen_vecs_true)
         N.testing.assert_allclose(correlation_mat, correlation_mat_true)
         N.testing.assert_allclose(modes, modes_true[:,self.mode_indices])
                     
         modes, eigen_vals, eigen_vecs = \
             compute_POD_matrices_direct_method(vec_array, self.mode_indices, 
             inner_product_weights=weights, return_all=True)
         
         N.testing.assert_allclose(eigen_vals, eigen_vals_true)
         N.testing.assert_allclose(N.abs(eigen_vecs), N.abs(eigen_vecs_true))
         N.testing.assert_allclose(N.abs(modes), N.abs(modes_true[:,self.mode_indices]))

    def prob(self, state):
        """
        Calculates the probability of a given state using the exponential
        family parameterization learned from a prior dataset. Note that
        we are assuming A(theta) = 0 (i.e., the graph is triangulated).
        Also note the similarity to calc_mu.
        """
        n1 = state[0,0:self.num_sites]
        n2 = state[0,self.num_sites:]
        mu_s = state
        s = np.zeros((self.num_nodes, self.num_nodes))
        mu_st11 = np.matrix(s)
        mu_st10 = np.matrix(s)
        mu_st01 = np.matrix(s)
        mu_st00 = np.matrix(s)

        """
        Calculate the edges from n->n'.
        This is the upper-right quadrant.
        """	
        upright = self.calc_mu_quadrant(n1, n2, False)
        mu_st11[0:self.num_sites,self.num_sites:] += upright[3] # nn11
        mu_st10[0:self.num_sites,self.num_sites:] += upright[2] # nn10
        mu_st01[0:self.num_sites,self.num_sites:] += upright[1] # nn01
        mu_st00[0:self.num_sites,self.num_sites:] += upright[0] # nn00
        """
        Calculate the edges from n'->n'.
        This is the lower-right quadrant.
        """
        lowright = self.calc_mu_quadrant(n2, n2, True)
        mu_st11[self.num_sites:,self.num_sites:] += lowright[3] # nn11
        mu_st10[self.num_sites:,self.num_sites:] += lowright[2] # nn10
        mu_st01[self.num_sites:,self.num_sites:] += lowright[1] # nn01
        mu_st00[self.num_sites:,self.num_sites:] += lowright[0] # nn00
        """
        We only want the upper triangle, since otherwise we will double
        count the n'->n' edges.
        """
        mu_st11 = np.triu(mu_st11)
        mu_st10 = np.triu(mu_st10)
        mu_st01 = np.triu(mu_st01)
        mu_st00 = np.triu(mu_st00)
        """
        Sum(x_s)
        """
        sum_s = mu_s * self.nodes[0].T
        sum_s += (mu_s - 1) * -1 * self.nodes[1].T
        sum_s = sum_s[0,0]
        """
        Sum(x_s,x_t)
        """
        sum_st00 = np.sum(np.multiply(mu_st00, self.edges[0]))
        sum_st01 = np.sum(np.multiply(mu_st01, self.edges[1]))
        sum_st10 = np.sum(np.multiply(mu_st10, self.edges[2]))
        sum_st11 = np.sum(np.multiply(mu_st11, self.edges[3]))
        result = sum_s + sum_st00 + sum_st01 + sum_st10 + sum_st11
        #print "S: {0} ST[00]: {1} ST[01]: {2} ST[10]: {3} ST[11]: {4}".format(sum_s, sum_st00, sum_st01, sum_st10, sum_st11)
        #print "Prob: {0}".format(result)
        #pseudo-likelihood
        return result

def make_connections(n,density=0.25):
	"""
	This function will return a random adjacency matrix of size
	n x n. You read the matrix like this:
	
	if matrix[2,7] = 1, then cities '2' and '7' are connected.
	if matrix[2,7] = 0, then the cities are _not_ connected.
	
	:param n: number of cities
	:param density: controls the ratio of 1s to 0s in the matrix
	
	:returns: an n x n adjacency matrix
	"""
	
	import networkx
	
	# Generate a random adjacency matrix and use it to build a networkx graph
	a=numpy.int32(numpy.triu((numpy.random.random_sample(size=(n,n))<density)))
	G=networkx.from_numpy_matrix(a)
	
	# If the network is 'not connected' (i.e., there are isolated nodes)
	# generate a new one. Keep doing this until we get a connected one.
	# Yes, there are more elegant ways to do this, but I'm demonstrating
	# while loops!
	while not networkx.is_connected(G):
		a=numpy.int32(numpy.triu((numpy.random.random_sample(size=(n,n))<density)))
		G=networkx.from_numpy_matrix(a)
	
	# Cities should be connected to themselves.
	numpy.fill_diagonal(a,1)
	
	return a + numpy.triu(a,1).T

    def __init__(self, rng, matches_vec, batch_size,
            sample_diff_every_epoch=True, n_same_pairs=None):
        """
        If `n_same_pairs` is given, this number of same pairs is sampled,
        otherwise all same pairs are used.
        """
        self.rng = rng
        self.matches_vec = matches_vec
        self.batch_size = batch_size

        if n_same_pairs is None:
            # Use all pairs
            I, J = np.where(np.triu(distance.squareform(matches_vec)))  # indices of same pairs
        else:
            # Sample same pairs
            n_pairs = min(n_same_pairs, len(np.where(matches_vec == True)[0]))
            same_sample = self.rng.choice(
                np.where(matches_vec == True)[0], size=n_pairs, replace=False
                )
            same_vec = np.zeros(self.matches_vec.shape[0], dtype=np.bool)
            same_vec[same_sample] = True
            I, J = np.where(np.triu(distance.squareform(same_vec)))

        self.x1_same_indices = []
        self.x2_same_indices = []
        for i, j in zip(I, J):
            self.x1_same_indices.append(i)
            self.x2_same_indices.append(j)

        if not sample_diff_every_epoch:
            self.x1_diff_indices, self.x2_diff_indices = self._sample_diff_pairs()
        self.sample_diff_every_epoch = sample_diff_every_epoch

    def __call__(self, container, sim_wide_params):
        distance_matrices = container.get_distance_matrices()
        # LJ
        accelerations, pe = moldyn.lennardJonesForce(
                distance_matrices,
                radius=self.radius,
                epsilon=self.epsilon,
                cutoff_dist=self.cutoff_dist,
                anchor_ixs = sim_wide_params.anchor_ixs)
        self.pe = pe  # should LJ func calc this differently since we are doing a cutoff?
        # Damping
        dx, dy, dr = distance_matrices  # KLUDGE
        ixs_not_touching = dr > self.cutoff_dist
        self.ixs_not_touching = ixs_not_touching
        dvx, dvy, dvr = moldyn.get_distance_matrices(container.velocities)
##        if np.any(container.velocities > 1):
##            pass
##        if container.time > 1:
##            pass
        dr[dr == 0] = 0.1  # prevent / 0 errors. Will this cause problems?
        accel = self.damping_magnitude * (dvx*dx + dvy*dy)/dr
        accel[ixs_not_touching] = 0
        accelerations[:, 0] += np.sum(np.triu(accel * dx/dr), axis=0)
        accelerations[:, 1] += np.sum(np.triu(accel * dy/dr), axis=0)
        return accelerations