Python numpy.diag_indices() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def test_covariances_and_eigenvalues(self):
        reader = FeatureReader(self.trajnames, self.temppdb)
        for tau in [1, 10, 100, 1000, 2000]:
            trans = tica(lag=tau, dim=self.dim, kinetic_map=False)
            trans.data_producer = reader

            log.info('number of trajectories reported by tica %d' % trans.number_of_trajectories())
            trans.parametrize()
            data = trans.get_output()

            log.info('max. eigenvalue: %f' % np.max(trans.eigenvalues))
            self.assertTrue(np.all(trans.eigenvalues <= 1.0))
            # check ICs
            check = tica(data=data, lag=tau, dim=self.dim)

            np.testing.assert_allclose(np.eye(self.dim), check.cov, atol=1e-8)
            np.testing.assert_allclose(check.mean, 0.0, atol=1e-8)
            ic_cov_tau = np.zeros((self.dim, self.dim))
            ic_cov_tau[np.diag_indices(self.dim)] = trans.eigenvalues
            np.testing.assert_allclose(ic_cov_tau, check.cov_tau, atol=1e-8) 

Example 2

def update_hypers(self, params):
        for i in range(self.nolayers):
            layeri = self.layers[i]
            Mi = layeri.M
            Dini = layeri.Din
            Douti = layeri.Dout
            layeri.ls.set_value(params['ls'][i])
            layeri.sf.set_value(params['sf'][i])
            layeri.sn.set_value(params['sn'][i])
            triu_ind = np.triu_indices(Mi)
            diag_ind = np.diag_indices(Mi)
            for d in range(Douti):
                layeri.zu[d].set_value(params['zu'][i][d])
                theta_m_d = params['eta2'][i][d]
                theta_R_d = params['eta1_R'][i][d]
                R = np.zeros((Mi, Mi))
                R[triu_ind] = theta_R_d.reshape(theta_R_d.shape[0], )
                R[diag_ind] = np.exp(R[diag_ind])
                layeri.theta_1_R[d] = R
                layeri.theta_1[d] = np.dot(R.T, R)
                layeri.theta_2[d] = theta_m_d 

Example 3

def _passthrough_delay(m, c):
    """Analytically derived state-space when p = q = m.

    We use this because it is numerically stable for high m.
    """
    j = np.arange(1, m+1, dtype=np.float64)
    u = (m + j) * (m - j + 1) / (c * j)

    A = np.zeros((m, m))
    B = np.zeros((m, 1))
    C = np.zeros((1, m))
    D = np.zeros((1,))

    A[0, :] = B[0, 0] = -u[0]
    A[1:, :-1][np.diag_indices(m-1)] = u[1:]
    D[0] = (-1)**m
    C[0, np.arange(m) % 2 == 0] = 2*D[0]
    return LinearSystem((A, B, C, D), analog=True) 

Example 4

def _proper_delay(q, c):
    """Analytically derived state-space when p = q - 1.

    We use this because it is numerically stable for high q
    and doesn't have a passthrough.
    """
    j = np.arange(q, dtype=np.float64)
    u = (q + j) * (q - j) / (c * (j + 1))

    A = np.zeros((q, q))
    B = np.zeros((q, 1))
    C = np.zeros((1, q))
    D = np.zeros((1,))

    A[0, :] = -u[0]
    B[0, 0] = u[0]
    A[1:, :-1][np.diag_indices(q-1)] = u[1:]
    C[0, :] = (j + 1) / float(q) * (-1) ** (q - 1 - j)
    return LinearSystem((A, B, C, D), analog=True) 

Example 5

def convert_solution(z, Cs):
    if issparse(Cs):
        Cs = Cs.toarray()
    M = Cs.shape[0]
    x = z[0:M]
    y = z[M:]

    w=np.exp(y)
    pi=w/w.sum()
    
    X=pi[:,np.newaxis]*x[np.newaxis,:]
    Y=X+np.transpose(X)
    denom=Y
    enum=Cs*np.transpose(pi)
    P=enum/denom
    ind=np.diag_indices(Cs.shape[0])
    P[ind]=0.0
    rowsums=P.sum(axis=1)
    P[ind]=1.0-rowsums
    return pi, P

###############################################################################
# Objective, Gradient, and Hessian
############################################################################### 

Example 6

def get_connectivity_matrix_nodiag(self):
        """
        Returns a similar matrix as in Ontology.get_connectivity_matrix(),
        but the diagonal of the matrix is 0.

        Note: !!!!!!!!!!!!!!!!!!!!!!!!
            d[a, a] == 0 instead of 1

        """

        if not hasattr(self, 'd_nodiag'):
            d = self.get_connectivity_matrix()
            self.d_nodiag = d.copy()
            self.d_nodiag[np.diag_indices(self.d_nodiag.shape[0])] = 0
            assert not np.isfortran(self.d_nodiag)

        return self.d_nodiag 

Example 7

def corrsort(features, use_tsp=False):
    """
    Given a features array, one feature per axis=0 entry, return axis=0 indices
    such that adjacent indices correspond to features that are correlated.

    cf. Traveling Salesman Problem. Not an optimal solution.

    use_tsp:
    Use tsp solver. See tsp_solver.greedy module that is used for this.
    Slows run-time considerably: O(N^4) computation, O(N^2) memory.

    Without use_tsp, both computation and memory are O(N^2).
    """

    correlations = np.ma.corrcoef(arrays.plane(features))
    if use_tsp: return tsp.solve_tsp(-correlations)

    size = features.shape[0]    
    correlations.mask[np.diag_indices(size)] = True
    
    # initialize results with the pair with the highest correlations.
    largest = np.argmax(correlations)
    results = [int(largest / size), largest % size]
    correlations.mask[:, results[0]] = True

    while len(results) < size:
        correlations.mask[:, results[-1]] = True
        results.append(np.argmax(correlations[results[-1]]))
            
    return results 

Example 8

def _get_batch_diagonal_cpu(array):
    batch_size, m, n = array.shape
    assert m == n
    rows, cols = np.diag_indices(n)
    return array[:, rows, cols] 

Example 9

def _set_batch_diagonal_cpu(array, diag_val):
    batch_size, m, n = array.shape
    assert m == n
    rows, cols = np.diag_indices(n)
    array[:, rows, cols] = diag_val 

Example 10

def _diagonal_idx_array(batch_size, n):
    idx_offsets = np.arange(
        start=0, stop=batch_size * n * n, step=n * n, dtype=np.int32).reshape(
        (batch_size, 1))
    idx = np.ravel_multi_index(
        np.diag_indices(n), (n, n)).reshape((1, n)).astype(np.int32)
    return cuda.to_gpu(idx + idx_offsets) 

Example 11

def check_forward(self, diag_data, non_diag_data):
        diag = chainer.Variable(diag_data)
        non_diag = chainer.Variable(non_diag_data)
        y = lower_triangular_matrix(diag, non_diag)

        correct_y = numpy.zeros(
            (self.batch_size, self.n, self.n), dtype=numpy.float32)

        tril_rows, tril_cols = numpy.tril_indices(self.n, -1)
        correct_y[:, tril_rows, tril_cols] = cuda.to_cpu(non_diag_data)

        diag_rows, diag_cols = numpy.diag_indices(self.n)
        correct_y[:, diag_rows, diag_cols] = cuda.to_cpu(diag_data)

        gradient_check.assert_allclose(correct_y, cuda.to_cpu(y.data)) 

Example 12

def kth_diag_indices(n, k):
    """
    Return indices of bins k steps away from the diagonal.
    (from http://stackoverflow.com/questions/10925671/numpy-k-th-diagonal-indices)
    """

    rows, cols = np.diag_indices(n)
    if k < 0:
        return rows[:k], cols[-k:]
    elif k > 0:
        return rows[k:], cols[:-k]
    else:
        return rows, cols 

Example 13

def compute_posterior_grad_u(self, dmu, dSu):
        # return grads wrt u params and Kuuinv
        triu_ind = np.triu_indices(self.M)
        diag_ind = np.diag_indices(self.M)
        if self.nat_param:
            dSu_via_m = np.einsum('da,db->dab', dmu, self.theta_2)
            dSu += dSu_via_m
            dSuinv = - np.einsum('dab,dbc,dce->dae', self.Su, dSu, self.Su)
            dKuuinv = np.sum(dSuinv, axis=0)
            dtheta1 = dSuinv
            deta2 = np.einsum('dab,db->da', self.Su, dmu)
        else:
            deta2 = dmu
            dtheta1 = dSu
            dKuuinv = 0

        dtheta1T = np.transpose(dtheta1, [0, 2, 1])
        dtheta1_R = np.einsum('dab,dbc->dac', self.theta_1_R, dtheta1 + dtheta1T)
        deta1_R = np.zeros([self.Dout, self.M * (self.M + 1) / 2])
        for d in range(self.Dout):
            dtheta1_R_d = dtheta1_R[d, :, :]
            theta1_R_d = self.theta_1_R[d, :, :]
            dtheta1_R_d[diag_ind] = dtheta1_R_d[diag_ind] * theta1_R_d[diag_ind]
            dtheta1_R_d = dtheta1_R_d[triu_ind]
            deta1_R[d, :] = dtheta1_R_d.reshape((dtheta1_R_d.shape[0], ))

        return deta1_R, deta2, dKuuinv 

Example 14

def get_hypers(self, key_suffix=''):
        """Summary
        
        Args:
            key_suffix (str, optional): Description
        
        Returns:
            TYPE: Description
        """
        params = {}
        M = self.M
        Din = self.Din
        Dout = self.Dout
        params['ls' + key_suffix] = self.ls
        params['sf' + key_suffix] = self.sf
        triu_ind = np.triu_indices(M)
        diag_ind = np.diag_indices(M)
        params_eta2 = self.theta_2
        params_eta1_R = np.zeros((Dout, M * (M + 1) / 2))
        params_zu_i = self.zu

        for d in range(Dout):
            Rd = np.copy(self.theta_1_R[d, :, :])
            Rd[diag_ind] = np.log(Rd[diag_ind])
            params_eta1_R[d, :] = np.copy(Rd[triu_ind])

        params['zu' + key_suffix] = self.zu
        params['eta1_R' + key_suffix] = params_eta1_R
        params['eta2' + key_suffix] = params_eta2
        return params 

Example 15

def update_hypers(self, params, key_suffix=''):
        """Summary
        
        Args:
            params (TYPE): Description
            key_suffix (str, optional): Description
        """
        M = self.M
        self.ls = params['ls' + key_suffix]
        self.sf = params['sf' + key_suffix]
        triu_ind = np.triu_indices(M)
        diag_ind = np.diag_indices(M)
        zu = params['zu' + key_suffix]
        self.zu = zu

        for d in range(self.Dout):
            theta_m_d = params['eta2' + key_suffix][d, :]
            theta_R_d = params['eta1_R' + key_suffix][d, :]
            R = np.zeros((M, M))
            R[triu_ind] = np.copy(theta_R_d.reshape(theta_R_d.shape[0], ))
            R[diag_ind] = np.exp(R[diag_ind])
            self.theta_1_R[d, :, :] = R
            self.theta_1[d, :, :] = np.dot(R.T, R)
            self.theta_2[d, :] = theta_m_d

        # update Kuu given new hypers
        self.compute_kuu()
        # compute mu and Su for each layer
        self.update_posterior() 

Example 16

def compute_cav_grad_u(self, dmu, dSu, alpha):
        # return grads wrt u params and Kuuinv
        triu_ind = np.triu_indices(self.M)
        diag_ind = np.diag_indices(self.M)
        beta = (self.N - alpha) * 1.0 / self.N
        if self.nat_param:
            dSu_via_m = np.einsum('da,db->dab', dmu, beta * self.theta_2)
            dSu += dSu_via_m
            dSuinv = - np.einsum('dab,dbc,dce->dae', self.Suhat, dSu, self.Suhat)
            dKuuinv = np.sum(dSuinv, axis=0)
            dtheta1 = beta * dSuinv
            deta2 = beta * np.einsum('dab,db->da', self.Suhat, dmu)
        else:
            Suhat = self.Suhat
            Suinv = self.Suinv
            mu = self.mu
            data_f_2 = np.einsum('dab,db->da', Suinv, mu)
            dSuhat_via_mhat = np.einsum('da,db->dab', dmu, beta * data_f_2)
            dSuhat = dSu + dSuhat_via_mhat
            dmuhat = dmu
            dSuhatinv = - np.einsum('dab,dbc,dce->dae', Suhat, dSuhat, Suhat)
            dSuinv_1 = beta * dSuhatinv
            Suhatdmu = np.einsum('dab,db->da', Suhat, dmuhat)
            dSuinv = dSuinv_1 + beta * np.einsum('da,db->dab', Suhatdmu, mu)
            dtheta1 = - np.einsum('dab,dbc,dce->dae', Suinv, dSuinv, Suinv)
            deta2 = beta * np.einsum('dab,db->da', Suinv, Suhatdmu)
            dKuuinv = (1 - beta) / beta * np.sum(dSuinv_1, axis=0)

        dtheta1T = np.transpose(dtheta1, [0, 2, 1])
        dtheta1_R = np.einsum('dab,dbc->dac', self.theta_1_R, dtheta1 + dtheta1T)
        deta1_R = np.zeros([self.Dout, self.M * (self.M + 1) / 2])
        for d in range(self.Dout):
            dtheta1_R_d = dtheta1_R[d, :, :]
            theta1_R_d = self.theta_1_R[d, :, :]
            dtheta1_R_d[diag_ind] = dtheta1_R_d[diag_ind] * theta1_R_d[diag_ind]
            dtheta1_R_d = dtheta1_R_d[triu_ind]
            deta1_R[d, :] = dtheta1_R_d.reshape((dtheta1_R_d.shape[0], ))

        return deta1_R, deta2, dKuuinv 

Example 17

def test_integration_adaptive_graph_lasso(self, params_in):
        '''
        Just tests inputs/outputs (not validity of result).
        '''
        n_features = 20
        n_samples = 25
        cov, prec, adj = ClusterGraph(
            n_blocks=1,
            chain_blocks=False,
            seed=1,
        ).create(n_features, 0.8)
        prng = np.random.RandomState(2)
        X = prng.multivariate_normal(np.zeros(n_features), cov, size=n_samples)

        model = AdaptiveGraphLasso(**params_in)
        model.fit(X)
        assert model.estimator_ is not None
        assert model.lam_ is not None

        assert np.sum(model.lam_[np.diag_indices(n_features)]) == 0

        if params_in['method'] == 'binary':
            uvals = set(model.lam_.flat)
            assert len(uvals) == 2
            assert 0 in uvals
            assert 1 in uvals
        elif params_in['method'] == 'inverse' or\
                params_in['method'] == 'inverse_squared':
            uvals = set(model.lam_.flat[model.lam_.flat != 0])
            assert len(uvals) > 0 

Example 18

def _nonzero_intersection(m, m_hat):
    '''Count the number of nonzeros in and between m and m_hat.

    Returns
    ----------
    m_nnz :  number of nonzeros in m (w/o diagonal)

    m_hat_nnz : number of nonzeros in m_hat (w/o diagonal)

    intersection_nnz : number of nonzeros in intersection of m/m_hat
                      (w/o diagonal)
    '''
    n_features, _ = m.shape

    m_no_diag = m.copy()
    m_no_diag[np.diag_indices(n_features)] = 0
    m_hat_no_diag = m_hat.copy()
    m_hat_no_diag[np.diag_indices(n_features)] = 0

    m_hat_nnz = len(np.nonzero(m_hat_no_diag.flat)[0])
    m_nnz = len(np.nonzero(m_no_diag.flat)[0])

    intersection_nnz = len(np.intersect1d(
        np.nonzero(m_no_diag.flat)[0],
        np.nonzero(m_hat_no_diag.flat)[0]
    ))

    return m_nnz, m_hat_nnz, intersection_nnz 

Example 19

def _binary_weights(self, estimator):
        n_features, _ = estimator.precision_.shape
        lam = np.zeros((n_features, n_features))
        lam[estimator.precision_ == 0] = 1
        lam[np.diag_indices(n_features)] = 0
        return lam 

Example 20

def _inverse_squared_weights(self, estimator):
        n_features, _ = estimator.precision_.shape
        lam = np.zeros((n_features, n_features))
        mask = estimator.precision_ != 0
        lam[mask] = 1. / (np.abs(estimator.precision_[mask]) ** 2)
        mask_0 = estimator.precision_ == 0
        lam[mask_0] = np.max(lam[mask].flat)  # non-zero in appropriate range
        lam[np.diag_indices(n_features)] = 0
        return lam 

Example 21

def _inverse_weights(self, estimator):
        n_features, _ = estimator.precision_.shape
        lam = np.zeros((n_features, n_features))
        mask = estimator.precision_ != 0
        lam[mask] = 1. / np.abs(estimator.precision_[mask])
        mask_0 = estimator.precision_ == 0
        lam[mask_0] = np.max(lam[mask].flat)  # non-zero in appropriate range
        lam[np.diag_indices(n_features)] = 0
        return lam 

Example 22

def _count_support_diff(m, m_hat):
    n_features, _ = m.shape

    m_no_diag = m.copy()
    m_no_diag[np.diag_indices(n_features)] = 0
    m_hat_no_diag = m_hat.copy()
    m_hat_no_diag[np.diag_indices(n_features)] = 0

    m_nnz = len(np.nonzero(m_no_diag.flat)[0])
    m_hat_nnz = len(np.nonzero(m_hat_no_diag.flat)[0])

    nnz_intersect = len(np.intersect1d(np.nonzero(m_no_diag.flat)[0],
                                       np.nonzero(m_hat_no_diag.flat)[0]))
    return m_nnz + m_hat_nnz - (2 * nnz_intersect) 

Example 23

def dist_diff(self, geom):
        diff = self.coords[:, None, :]-geom.coords[None, :, :]
        dist = np.sqrt(np.sum(diff**2, 2))
        dist[np.diag_indices(len(self))] = np.inf
        return dist, diff 

Example 24

def forward(self, x):
        # create diagonal matrices
        m = np.zeros((x.size * self.dim)).reshape(-1, self.dim, self.dim)
        x = x.reshape(-1, self.dim)
        m[(np.s_[:],) + np.diag_indices(x.shape[1])] = x
        return m 

Example 25

def __init__(self, n=10, l=None):

		self.n = n
		self.l = l

		self.W = randfilt(self.W if hasattr(self, 'W') else None, (self.l, self.n))

		self.C = np.zeros((n, n), dtype='float32')
		self.C[np.diag_indices(n)] = 1 

Example 26

def __init__(self, n=10, l=None):

		self.n = n
		self.l = l

		self.W = randfilt(self.W if hasattr(self, 'W') else None, (self.l, self.n))

		self.C = np.zeros((n, n), dtype='float32') - 1
		self.C[np.diag_indices(n)] = 1

		self.SII = None
		self.SIO = None 

Example 27

def test_naf_layer_full():
    batch_size = 2
    for nb_actions in (1, 3):
        # Construct single model with NAF as the only layer, hence it is fully deterministic
        # since no weights are used, which would be randomly initialized.
        L_flat_input = Input(shape=((nb_actions * nb_actions + nb_actions) // 2,))
        mu_input = Input(shape=(nb_actions,))
        action_input = Input(shape=(nb_actions,))
        x = NAFLayer(nb_actions, mode='full')([L_flat_input, mu_input, action_input])
        model = Model(inputs=[L_flat_input, mu_input, action_input], outputs=x)
        model.compile(loss='mse', optimizer='sgd')
        
        # Create random test data.
        L_flat = np.random.random((batch_size, (nb_actions * nb_actions + nb_actions) // 2)).astype('float32')
        mu = np.random.random((batch_size, nb_actions)).astype('float32')
        action = np.random.random((batch_size, nb_actions)).astype('float32')

        # Perform reference computations in numpy since these are much easier to verify.
        L = np.zeros((batch_size, nb_actions, nb_actions)).astype('float32')
        LT = np.copy(L)
        for l, l_T, l_flat in zip(L, LT, L_flat):
            l[np.tril_indices(nb_actions)] = l_flat
            l[np.diag_indices(nb_actions)] = np.exp(l[np.diag_indices(nb_actions)])
            l_T[:, :] = l.T
        P = np.array([np.dot(l, l_T) for l, l_T in zip(L, LT)]).astype('float32')
        A_ref = np.array([np.dot(np.dot(a - m, p), a - m) for a, m, p in zip(action, mu, P)]).astype('float32')
        A_ref *= -.5

        # Finally, compute the output of the net, which should be identical to the previously
        # computed reference.
        A_net = model.predict([L_flat, mu, action]).flatten()
        assert_allclose(A_net, A_ref, rtol=1e-5) 

Example 28

def test_naf_layer_diag():
    batch_size = 2
    for nb_actions in (1, 3):
        # Construct single model with NAF as the only layer, hence it is fully deterministic
        # since no weights are used, which would be randomly initialized.
        L_flat_input = Input(shape=(nb_actions,))
        mu_input = Input(shape=(nb_actions,))
        action_input = Input(shape=(nb_actions,))
        x = NAFLayer(nb_actions, mode='diag')([L_flat_input, mu_input, action_input])
        model = Model(inputs=[L_flat_input, mu_input, action_input], outputs=x)
        model.compile(loss='mse', optimizer='sgd')
        
        # Create random test data.
        L_flat = np.random.random((batch_size, nb_actions)).astype('float32')
        mu = np.random.random((batch_size, nb_actions)).astype('float32')
        action = np.random.random((batch_size, nb_actions)).astype('float32')

        # Perform reference computations in numpy since these are much easier to verify.
        P = np.zeros((batch_size, nb_actions, nb_actions)).astype('float32')
        for p, l_flat in zip(P, L_flat):
            p[np.diag_indices(nb_actions)] = l_flat
        print(P, L_flat)
        A_ref = np.array([np.dot(np.dot(a - m, p), a - m) for a, m, p in zip(action, mu, P)]).astype('float32')
        A_ref *= -.5

        # Finally, compute the output of the net, which should be identical to the previously
        # computed reference.
        A_net = model.predict([L_flat, mu, action]).flatten()
        assert_allclose(A_net, A_ref, rtol=1e-5) 

Example 29

def setupMatrices(dim, mult=0.05):
        di  = np.diag_indices(dim)
        bt  = np.random.randn(dim,)*mult
        Wt  = np.zeros((dim,dim))
        Wt[di] = np.random.randn(dim,)*mult
        We  = np.zeros((dim,dim))
        We[di] = np.random.randn(dim,)*mult
        return Wt,bt,We 

Example 30

def toPartialCorr(Prec):
    D=Prec[np.diag_indices(Prec.shape[0])]
    P=Prec.copy()
    D=np.outer(D,D)
    return -P/np.sqrt(D) 

Example 31

def toPartialCorr(Prec):
    D=Prec[np.diag_indices(Prec.shape[0])]
    P=Prec.copy()
    D=np.outer(D,D)
    return -P/np.sqrt(D) 

Example 32

def toPartialCorr(Prec):
    D=Prec[np.diag_indices(Prec.shape[0])]
    P=Prec.copy()
    D=np.outer(D,D)
    D[D == 0] = 1
    if(np.sum(D==0)):
        print 'Ds',np.sum(D==0)
    if(np.sum(D<0)):
        print 'Dsminus',np.sum(D<0)        
    return -P/np.sqrt(D) 

Example 33

def corrcoef_matrix(matrix):
    # Code originating from http://stackoverflow.com/a/24547964 by http://stackoverflow.com/users/2455058/jingchao

    r = np.corrcoef(matrix)
    rf = r[np.triu_indices(r.shape[0], 1)]
    df = matrix.shape[1] - 2
    ts = rf * rf * (df / (1 - rf * rf))
    pf = betai(0.5 * df, 0.5, df / (df + ts))
    p = np.zeros(shape=r.shape)
    p[np.triu_indices(p.shape[0], 1)] = pf
    p[np.tril_indices(p.shape[0], -1)] = pf
    p[np.diag_indices(p.shape[0])] = np.ones(p.shape[0])
    return r, p 

Example 34

def get_hypers(self):
        params = {'ls': [],
                  'sf': [],
                  'zu': [],
                  'sn': [],
                  'eta1_R': [],
                  'eta2': []}
        for i in range(self.nolayers):
            layeri = self.layers[i]
            Mi = layeri.M
            Dini = layeri.Din
            Douti = layeri.Dout
            params['ls'].append(layeri.ls.get_value())
            params['sf'].append(layeri.sf.get_value())
            params['sn'].append(layeri.sn.get_value())
            triu_ind = np.triu_indices(Mi)
            diag_ind = np.diag_indices(Mi)
            params_zu_i = []
            params_eta2_i = []
            params_eta1_Ri = []
            for d in range(Douti):
                params_zu_i.append(layeri.zu[d].get_value())
                params_eta2_i.append(layeri.theta_2[d])

                Rd = layeri.theta_1_R[d]
                Rd[diag_ind] = np.log(Rd[diag_ind])
                params_eta1_Ri.append(Rd[triu_ind].reshape((Mi*(Mi+1)/2, 1)))

            params['zu'].append(params_zu_i)
            params['eta1_R'].append(params_eta1_Ri)
            params['eta2'].append(params_eta2_i)
        return params 

Example 35

def get_hypers(self):
        params = {'ls': [],
                  'sf': [],
                  'zu': [],
                  'sn': [],
                  'eta1_R': [],
                  'eta2': []}
        for i in range(self.no_layers):
            Mi = self.no_pseudos[i]
            Dini = self.layer_sizes[i]
            Douti = self.layer_sizes[i+1]
            params['ls'].append(self.ls[i])
            params['sf'].append(self.sf[i])
            if not (self.no_output_noise and (i == self.no_layers-1)):
                params['sn'].append(self.sn[i])
            triu_ind = np.triu_indices(Mi)
            diag_ind = np.diag_indices(Mi)
            params_zu_i = []
            params_eta2_i = []
            params_eta1_Ri = []
            if self.zu_tied:
                params_zu_i = self.zu[i]
            else:
                for d in range(Douti):
                    params_zu_i.append(self.zu[i][d])

            for d in range(Douti):
                params_eta2_i.append(self.theta_2[i][d])
                Rd = self.theta_1_R[i][d]
                Rd[diag_ind] = np.log(Rd[diag_ind])
                params_eta1_Ri.append(Rd[triu_ind].reshape((Mi*(Mi+1)/2,)))

            params['zu'].append(params_zu_i)
            params['eta1_R'].append(params_eta1_Ri)
            params['eta2'].append(params_eta2_i)
        return params 

Example 36

def update_hypers(self, params):
        for i in range(self.no_layers):
            Mi = self.no_pseudos[i]
            Dini = self.layer_sizes[i]
            Douti = self.layer_sizes[i+1]
            self.ls[i] = params['ls'][i]
            self.ones_M_ls[i] = np.outer(self.ones_M[i], 1.0 / np.exp(2*self.ls[i]))
            self.sf[i] = params['sf'][i]
            if not ((self.no_output_noise) and (i == self.no_layers-1)):
                self.sn[i] = params['sn'][i]
            triu_ind = np.triu_indices(Mi)
            diag_ind = np.diag_indices(Mi)
            if self.zu_tied:
                zi = params['zu'][i]
                self.zu[i] = zi
            else:
                for d in range(Douti):
                    zid = params['zu'][i][d]
                    self.zu[i][d] = zid

            for d in range(Douti):
                theta_m_d = params['eta2'][i][d]
                theta_R_d = params['eta1_R'][i][d]
                R = np.zeros((Mi, Mi))
                R[triu_ind] = theta_R_d.reshape(theta_R_d.shape[0], )
                R[diag_ind] = np.exp(R[diag_ind])
                self.theta_1_R[i][d] = R
                self.theta_1[i][d] = np.dot(R.T, R)
                self.theta_2[i][d] = theta_m_d 

Example 37

def test_balreal():
    isys = Lowpass(0.05)
    noise = 0.5*Lowpass(0.01) + 0.5*Alpha(0.005)
    p = 0.8
    sys = p*isys + (1-p)*noise

    T, Tinv, S = balanced_transformation(sys)
    assert np.allclose(inv(T), Tinv)
    assert np.allclose(S, hsvd(sys))

    balsys = sys.transform(T, Tinv)
    assert balsys == sys

    assert np.all(S >= 0)
    assert np.all(S[0] > 0.3)
    assert np.all(S[1:] < 0.05)
    assert np.allclose(sorted(S, reverse=True), S)

    P = control_gram(balsys)
    Q = observe_gram(balsys)

    diag = np.diag_indices(len(P))
    offdiag = np.ones_like(P, dtype=bool)
    offdiag[diag] = False
    offdiag = np.where(offdiag)

    assert np.allclose(P[diag], S)
    assert np.allclose(P[offdiag], 0)
    assert np.allclose(Q[diag], S)
    assert np.allclose(Q[offdiag], 0) 

Example 38

def __add_third_in_line(matrix, player, win_or_defend):
    """Look for 2 in line and try to win or defend in 1 move"""
    done = False
    opponent = model.opponent(player)
    # = horizontals
    for row in range(3):
        if not done:
            done, line = win_or_defend(matrix[row, :], player)
            if done:
                matrix[row, :] = line
                # || verticals
    for col in range(3):
        if not done:
            done, line = win_or_defend(matrix[:, col], player)
            if done:
                matrix[:, col] = line
    if not done:
        # \ diagonal
        done, line = win_or_defend(matrix.diagonal(), player)
        if done:
            matrix[np.diag_indices(3)] = line
    if not done:
        # / diagonal
        done, line = win_or_defend(np.fliplr(matrix).diagonal(), player)
        if done:
            np.fliplr(matrix)[np.diag_indices(3)] = line
    if not done:
        return matrix, player
    else:
        return matrix, opponent 

Example 39

def _get_diagIDs(K):
  if K in diagIDsDict:
    return diagIDsDict[K]
  else:
    diagIDs = np.diag_indices(K)
    diagIDsDict[K] = diagIDs
    return diagIDs 

Example 40

def _calcWordTypeCooccurMatrix(self, Q, sameWordVec, nDoc):
        """ Transform building blocks into the final Q matrix

        Returns
        -------
        Q : 2D array, size W x W (where W is vocab_size)
        """
        Q /= np.float32(nDoc)
        sameWordVec /= np.float32(nDoc)
        diagIDs = np.diag_indices(self.vocab_size)
        Q[diagIDs] -= sameWordVec

        # Fix small numerical issues (like diag entries of -1e-15 instead of 0)
        np.maximum(Q, 1e-100, out=Q)
        return Q 

Example 41

def _calcWordTypeCooccurMatrix(self, Q, sameWordVec, nDoc):
        """ Transform building blocks into the final Q matrix

        Returns
        -------
        Q : 2D array, size W x W (where W is vocab_size)
        """
        Q /= np.float32(nDoc)
        sameWordVec /= np.float32(nDoc)
        diagIDs = np.diag_indices(self.vocab_size)
        Q[diagIDs] -= sameWordVec

        # Fix small numerical issues (like diag entries of -1e-15 instead of 0)
        np.maximum(Q, 1e-100, out=Q)
        return Q 

Example 42

def _get_diagIDs(K):
    if K in diagIDsDict:
        return diagIDsDict[K]
    else:
        diagIDs = np.diag_indices(K)
        diagIDsDict[K] = diagIDs
        return diagIDs 

Example 43

def _get_diagIDs(K):
    if K in diagIDsDict:
        return diagIDsDict[K]
    else:
        diagIDs = np.diag_indices(K)
        diagIDsDict[K] = diagIDs
        return diagIDs 

Example 44

def add_preference(hdf5_file, preference):
    """Assign the value 'preference' to the diagonal entries
        of the matrix of similarities stored in the HDF5 data structure 
        at 'hdf5_file'.
    """

    Worker.hdf5_lock.acquire()
    
    with tables.open_file(hdf5_file, 'r+') as fileh:
        S = fileh.root.aff_prop_group.similarities
        diag_ind = np.diag_indices(S.nrows)
        S[diag_ind] = preference
        
    Worker.hdf5_lock.release() 

Example 45

def check_convergence(hdf5_file, iteration, convergence_iter, max_iter):
    """If the estimated number of clusters has not changed for 'convergence_iter'
        consecutive iterations in a total of 'max_iter' rounds of message-passing, 
        the procedure herewith returns 'True'.
        Otherwise, returns 'False'.
        Parameter 'iteration' identifies the run of message-passing 
        that has just completed.
    """

    Worker.hdf5_lock.acquire()
    
    with tables.open_file(hdf5_file, 'r+') as fileh:
        A = fileh.root.aff_prop_group.availabilities
        R = fileh.root.aff_prop_group.responsibilities
        P = fileh.root.aff_prop_group.parallel_updates
        
        N = A.nrows
        diag_ind = np.diag_indices(N)
        
        E = (A[diag_ind] + R[diag_ind]) > 0
        P[:, iteration % convergence_iter] = E
        
        e_mat = P[:]
        K = E.sum(axis = 0)
        
    Worker.hdf5_lock.release()
        
    if iteration >= convergence_iter:
        se = e_mat.sum(axis = 1)
        unconverged = (np.sum((se == convergence_iter) + (se == 0)) != N)
        if (not unconverged and (K > 0)) or (iteration == max_iter):
                return True
                
    return False 

Example 46

def _gt_weights(W):
    """Computes the weights V for a Guttman transform V X = B(X) Z."""
    V = -W
    V[np.diag_indices(V.shape[0])] = W.sum(axis=1) - W.diagonal()

    return V 

Example 47

def _gt_mapping(D, W, Z):
    """Computes the mapping B(X) for a Guttman transform V X = B(X) Z."""
    # Compute the Euclidean distances between all pairs of points
    Dz = distance.cdist(Z, Z)
    # Fill the diagonal of Dz, because *we don't want a division by zero*
    np.fill_diagonal(Dz, 1e-5)

    B = - W * D / Dz
    np.fill_diagonal(B, 0.0)
    B[np.diag_indices(B.shape[0])] = -np.sum(B, axis=1)

    return B 

Example 48

def binning_as_matrix(
        bin_count,
        array,
        minimum=None,
        maximum=None,
        axis=[],
        remove_small_cycles=True):
    """
    :param bin_count:
    :param array:
    :param maximum: maximum value to be recognized. Values bigger than max
    will be filtered
    :param minimum: minimum value to be recognized. Values smaller than min
    will be filtered
    :param axis: list of placements for axis. Possible list-elements are:
            * 'bottom'
            * 'left'
            * 'right'
            * 'top'
    :param remove_small_cycles: if True cycles where start and end are
        identical after binning will be removed
    :return: data matrix with start in rows and target in columns
    """

    # Bilding a 2d-histogram

    if minimum is None:
        minimum = np.nanmin(array)
    if maximum is None:
        maximum = np.nanmax(array)

    minimum = float(minimum)
    maximum = float(maximum)

    start = array[:, 0]
    target = array[:, 1]

    classified_data = np.histogram2d(start, target, bins=bin_count, range=[
                                     [minimum, maximum], [minimum, maximum]])

    res_matrix = classified_data[0]
    intervall_edges = classified_data[1]

    axis_value = np.diff(intervall_edges) / 2.0 + intervall_edges[0:-1]

    if remove_small_cycles:
        indexes = np.diag_indices(bin_count)
        res_matrix[indexes] = 0.0

    if axis:
        res_matrix = append_axis(res_matrix, horizontal_axis=axis_value,
                                 vertical_axis=axis_value, axis=axis)

    return res_matrix 

Example 49

def get_displacement_tensor(self):
        """A matrix where the entry A[i, j, :] is the vector
        self.cartesian_pos[i] - self.cartesian_pos[j].

        For periodic systems the distance of an atom from itself is the
        smallest displacement of an atom from one of it's periodic copies, and
        the distance of two different atoms is the distance of two closest
        copies.

        Returns:
            np.array: 3D matrix containing the pairwise distance vectors.
        """
        if self._displacement_tensor is None:

            if self.pbc.any():
                pos = self.get_scaled_positions()
                disp_tensor = pos[:, None, :] - pos[None, :, :]

                # Take periodicity into account by wrapping coordinate elements
                # that are bigger than 0.5 or smaller than -0.5
                indices = np.where(disp_tensor > 0.5)
                disp_tensor[indices] = 1 - disp_tensor[indices]
                indices = np.where(disp_tensor < -0.5)
                disp_tensor[indices] = disp_tensor[indices] + 1

                # Transform to cartesian
                disp_tensor = self.to_cartesian(disp_tensor)

                # Figure out the smallest basis vector and set it as
                # displacement for diagonal
                cell = self.get_cell()
                basis_lengths = np.linalg.norm(cell, axis=1)
                min_index = np.argmin(basis_lengths)
                min_basis = cell[min_index]
                diag_indices = np.diag_indices(len(disp_tensor))
                disp_tensor[diag_indices] = min_basis

            else:
                pos = self.get_positions()
                disp_tensor = pos[:, None, :] - pos[None, :, :]

            self._displacement_tensor = disp_tensor

        return self._displacement_tensor 

Example 50

def test_parameterized_orf(self):
        T1 = 3.16e8
        pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12),
                            gamma=parameter.Uniform(1,7))
        orf = hd_orf_generic(a=parameter.Uniform(0,5),
                             b=parameter.Uniform(0,5),
                             c=parameter.Uniform(0,5))
        rn = gp_signals.FourierBasisGP(spectrum=pl, Tspan=T1, components=30)
        crn = gp_signals.FourierBasisCommonGP(spectrum=pl, orf=orf,
                                              components=30, name='gw',
                                              Tspan=T1)

        model = rn + crn
        pta = model(self.psrs[0]) + model(self.psrs[1])

        lA1, gamma1 = -13, 1e-15
        lA2, gamma2 = -13.3, 1e-15
        lAc, gammac = -13.1, 1e-15
        a, b, c = 1.9, 0.4, 0.23

        params = {'gw_log10_A': lAc, 'gw_gamma': gammac,
                  'gw_a': a, 'gw_b':b, 'gw_c':c,
                  'B1855+09_log10_A': lA1, 'B1855+09_gamma': gamma1,
                  'J1909-3744_log10_A': lA2, 'J1909-3744_gamma': gamma2}

        phi = pta.get_phi(params)
        phiinv = pta.get_phiinv(params)

        F1, f1 = utils.createfourierdesignmatrix_red(
            self.psrs[0].toas, nmodes=30, Tspan=T1)
        F2, f2 = utils.createfourierdesignmatrix_red(
            self.psrs[1].toas, nmodes=30, Tspan=T1)

        msg = 'F matrix incorrect'
        assert np.allclose(pta.get_basis(params)[0], F1, rtol=1e-10), msg
        assert np.allclose(pta.get_basis(params)[1], F2, rtol=1e-10), msg

        nftot = 120
        phidiag = np.zeros(nftot)
        phit = np.zeros((nftot, nftot))

        phidiag[:60] = utils.powerlaw(f1, lA1, gamma1)
        phidiag[:60] += utils.powerlaw(f1, lAc, gammac)
        phidiag[60:] = utils.powerlaw(f2, lA2, gamma2)
        phidiag[60:] += utils.powerlaw(f2, lAc, gammac)

        phit[np.diag_indices(nftot)] = phidiag
        orf = hd_orf_generic(self.psrs[0].pos, self.psrs[1].pos, a=a, b=b, c=c)
        spec = utils.powerlaw(f1, log10_A=lAc, gamma=gammac)
        phit[:60, 60:] = np.diag(orf*spec)
        phit[60:, :60] = phit[:60, 60:]

        msg = '{} {}'.format(np.diag(phi), np.diag(phit))
        assert np.allclose(phi, phit, rtol=1e-15, atol=1e-17), msg
        msg = 'PTA Phi inverse is incorrect {}.'.format(params)
        assert np.allclose(phiinv, np.linalg.inv(phit),
                           rtol=1e-15, atol=1e-17), msg 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注