Python numpy.diag_indices_from() Examples

The following are 30 code examples for showing how to use numpy.diag_indices_from(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: vnpy_crypto   Author: birforce   File: correlation_structures.py    License: MIT License 6 votes vote down vote up
def corr_equi(k_vars, rho):
    '''create equicorrelated correlation matrix with rho on off diagonal

    Parameters
    ----------
    k_vars : int
        number of variables, correlation matrix will be (k_vars, k_vars)
    rho : float
        correlation between any two random variables

    Returns
    -------
    corr : ndarray (k_vars, k_vars)
        correlation matrix

    '''
    corr = np.empty((k_vars, k_vars))
    corr.fill(rho)
    corr[np.diag_indices_from(corr)] = 1
    return corr 
Example 2
Project: pyGeoStatistics   Author: whimian   File: krige3d.py    License: MIT License 6 votes vote down vote up
def block_covariance(self):
        "return average covariance within block"
        if self._block_covariance is None:
            if self.ndb <= 1:  # point kriging
                self._block_covariance = self.unbias
            else:
                cov = list()
                for x1, y1, z1 in zip(self.xdb, self.ydb, self.zdb):
                    for x2, y2, z2 in zip(self.xdb, self.ydb, self.zdb):
                        # cov.append(self._cova3((x1, y1, z1), (x2, y2, z2)))
                        cov.append(cova3(
                            (x1, y1, z1), (x2, y2, z2),
                            self.rotmat, self.maxcov, self.nst,
                            self.it, self.cc, self.aa_hmax))
                cov = np.array(cov).reshape((self.ndb, self.ndb))
                cov[np.diag_indices_from(cov)] -= self.c0
                self._block_covariance = np.mean(cov)
        return self._block_covariance 
Example 3
Project: cupy   Author: cupy   File: insert.py    License: MIT License 6 votes vote down vote up
def diag_indices_from(arr):
    """
    Return the indices to access the main diagonal of an n-dimensional array.
    See `diag_indices` for full details.

    Args:
        arr (cupy.ndarray): At least 2-D.

     .. seealso:: :func:`numpy.diag_indices_from`

    """
    if not isinstance(arr, cupy.ndarray):
        raise TypeError("Argument must be cupy.ndarray")

    if not arr.ndim >= 2:
        raise ValueError("input array must be at least 2-d")
    # For more than d=2, the strided formula is only valid for arrays with
    # all dimensions equal, so we check first.
    if not cupy.all(cupy.diff(arr.shape) == 0):
        raise ValueError("All dimensions of input must be of equal length")

    return diag_indices(arr.shape[0], arr.ndim) 
Example 4
Project: nni   Author: microsoft   File: bayesian.py    License: MIT License 6 votes vote down vote up
def first_fit(self, train_x, train_y):
        """ Fit the regressor for the first time. """
        train_x, train_y = np.array(train_x), np.array(train_y)

        self._x = np.copy(train_x)
        self._y = np.copy(train_y)

        self._distance_matrix = edit_distance_matrix(self._x)
        k_matrix = bourgain_embedding_matrix(self._distance_matrix)
        k_matrix[np.diag_indices_from(k_matrix)] += self.alpha

        self._l_matrix = cholesky(k_matrix, lower=True)  # Line 2

        self._alpha_vector = cho_solve(
            (self._l_matrix, True), self._y)  # Line 3

        self._first_fitted = True
        return self 
Example 5
Project: q2mm   Author: ericchansen   File: gradient.py    License: MIT License 6 votes vote down vote up
def do_lagrange(ma, vb, factor):
    """
    Lagrange multipliers.

    Parameters
    ----------
    ma : NumPy matrix
    vb : NumPy vector
    factor : float
    """
    mac = copy.deepcopy(ma)
    ind = np.diag_indices_from(mac)
    logger.log(5, 'A:\n{}'.format(mac))
    mac[ind] = mac[ind] + factor
    logger.log(5, 'A:\n{}'.format(mac))
    changes = solver(mac, vb)
    return [('LAGRANGE F{}'.format(factor), changes)] 
Example 6
Project: q2mm   Author: ericchansen   File: gradient.py    License: MIT License 6 votes vote down vote up
def do_levenberg(ma, vb, factor):
    """
    Lagrange multipliers.

    Parameters
    ----------
    ma : NumPy matrix
    vb : NumPy vector
    factor : float
    """
    mac = copy.deepcopy(ma)
    ind = np.diag_indices_from(mac)
    mac[ind] = mac[ind] + factor
    logger.log(5, 'A:\n{}'.format(mac))
    changes = solver(mac, vb)
    return [('LM {}'.format(factor), changes)] 
Example 7
def corr_equi(k_vars, rho):
    '''create equicorrelated correlation matrix with rho on off diagonal

    Parameters
    ----------
    k_vars : int
        number of variables, correlation matrix will be (k_vars, k_vars)
    rho : float
        correlation between any two random variables

    Returns
    -------
    corr : ndarray (k_vars, k_vars)
        correlation matrix

    '''
    corr = np.empty((k_vars, k_vars))
    corr.fill(rho)
    corr[np.diag_indices_from(corr)] = 1
    return corr 
Example 8
Project: pyins   Author: nmayorov   File: filt.py    License: MIT License 6 votes vote down vote up
def _kalman_correct(x, P, z, H, R, gain_factor, gain_curve):
    PHT = np.dot(P, H.T)

    S = np.dot(H, PHT) + R
    e = z - H.dot(x)
    L = cholesky(S, lower=True)
    inn = solve_triangular(L, e, lower=True)

    if gain_curve is not None:
        q = (np.dot(inn, inn) / inn.shape[0]) ** 0.5
        f = gain_curve(q)
        if f == 0:
            return inn
        L *= (q / f) ** 0.5

    K = cho_solve((L, True), PHT.T, overwrite_b=True).T
    if gain_factor is not None:
        K *= gain_factor[:, None]

    U = -K.dot(H)
    U[np.diag_indices_from(U)] += 1
    x += K.dot(z - H.dot(x))
    P[:] = U.dot(P).dot(U.T) + K.dot(R).dot(K.T)

    return inn 
Example 9
Project: meanderpy   Author: zsylvester   File: meanderpy.py    License: Apache License 2.0 6 votes vote down vote up
def order_cl_pixels(x_pix,y_pix):
    dist = distance.cdist(np.array([x_pix,y_pix]).T,np.array([x_pix,y_pix]).T)
    dist[np.diag_indices_from(dist)]=100.0
    ind = np.argmin(x_pix) # select starting point on left side of image
    clinds = [ind]
    count = 0
    while count<len(x_pix):
        t = dist[ind,:].copy()
        if len(clinds)>2:
            t[clinds[-2]]=t[clinds[-2]]+100.0
            t[clinds[-3]]=t[clinds[-3]]+100.0
        ind = np.argmin(t)
        clinds.append(ind)
        count=count+1
    x_pix = x_pix[clinds]
    y_pix = y_pix[clinds]
    return x_pix,y_pix 
Example 10
Project: enspara   Author: bowman-lab   File: test_mutual_info.py    License: GNU General Public License v3.0 6 votes vote down vote up
def test_symmetrical_mi_nonzero():
    # test that the MI matrix for sets of uncorrelated things results
    # in zero MI

    nonzero_mi_funcs = [nonzero_mi_np, nonzero_mi_ra, nonzero_mi_list]
    for a, n_states in (f() for f in nonzero_mi_funcs):

        mi = mutual_info.mi_matrix(a, a, n_states, n_states)

        assert_almost_equal(mi[-1, -2], 0.86114, decimal=3)
        mi[-1, -2] = mi[-2, -1] = 0

        assert_almost_equal(np.diag(mi), 0.86114, decimal=2)
        mi[np.diag_indices_from(mi)] = 0

        assert_allclose(mi, 0, atol=1e-3) 
Example 11
Project: enspara   Author: bowman-lab   File: test_mutual_info.py    License: GNU General Public License v3.0 6 votes vote down vote up
def test_symmetrical_mi_nonzero_int_shape_spec():
    # test that when we use an integer (rather than a list of integers)
    # that we correctly assume that the integer is just repeated for all
    # of the various features.

    nonzero_mi_funcs = [nonzero_mi_np, nonzero_mi_ra, nonzero_mi_list]
    for a, n_states in (f() for f in nonzero_mi_funcs):

        mi = mutual_info.mi_matrix(a, a, 5, 5)

        assert_almost_equal(mi[-1, -2], 0.86114, decimal=3)
        mi[-1, -2] = mi[-2, -1] = 0

        assert_almost_equal(np.diag(mi), 0.86114, decimal=2)
        mi[np.diag_indices_from(mi)] = 0

        assert_allclose(mi, 0, atol=1e-3) 
Example 12
Project: pyscf   Author: pyscf   File: mole.py    License: Apache License 2.0 5 votes vote down vote up
def energy_nuc(mol, charges=None, coords=None):
    '''Compute nuclear repulsion energy (AU) or static Coulomb energy

    Returns
        float
    '''
    if charges is None: charges = mol.atom_charges()
    if len(charges) <= 1:
        return 0
    #e = 0
    #for j in range(len(mol._atm)):
    #    q2 = charges[j]
    #    r2 = coords[j]
    #    for i in range(j):
    #        q1 = charges[i]
    #        r1 = coords[i]
    #        r = numpy.linalg.norm(r1-r2)
    #        e += q1 * q2 / r
    rr = inter_distance(mol, coords)
    rr[numpy.diag_indices_from(rr)] = 1e200
    if CHECK_GEOM and numpy.any(rr < 1e-5):
        for atm_idx in numpy.argwhere(rr<1e-5):
            logger.warn(mol, 'Atoms %s have the same coordinates', atm_idx)
        raise RuntimeError('Ill geometry')
    e = numpy.einsum('i,ij,j->', charges, 1./rr, charges) * .5
    return e 
Example 13
Project: pyscf   Author: pyscf   File: mole.py    License: Apache License 2.0 5 votes vote down vote up
def inter_distance(mol, coords=None):
    '''
    Inter-particle distance array
    '''
    if coords is None: coords = mol.atom_coords()
    rr = numpy.linalg.norm(coords.reshape(-1,1,3) - coords, axis=2)
    rr[numpy.diag_indices_from(rr)] = 0
    return rr 
Example 14
Project: pyGSTi   Author: pyGSTio   File: matrixtools.py    License: Apache License 2.0 5 votes vote down vote up
def matrix_sign(M):
    """ The "sign" matrix of `M` """
    #Notes: sign(M) defined s.t. eigvecs of sign(M) are evecs of M
    # and evals of sign(M) are +/-1 or 0 based on sign of eigenvalues of M

    #Using the extremely numerically stable (but expensive) Schur method
    # see http://www.maths.manchester.ac.uk/~higham/fm/OT104HighamChapter5.pdf
    N = M.shape[0]; assert(M.shape == (N, N)), "M must be square!"
    T, Z = _spl.schur(M, 'complex')  # M = Z T Z^H where Z is unitary and T is upper-triangular
    U = _np.zeros(T.shape, 'complex')  # will be sign(T), which is easy to compute
    # (U is also upper triangular), and then sign(M) = Z U Z^H

    # diagonals are easy
    U[_np.diag_indices_from(U)] = _np.sign(_np.diagonal(T))

    #Off diagonals: use U^2 = I or TU = UT
    # Note: Tij = Uij = 0 when i > j and i==j easy so just consider i<j case
    # 0 = sum_k Uik Ukj =  (i!=j b/c off-diag)
    # FUTURE: speed this up by using np.dot instead of sums below
    for j in range(1, N):
        for i in range(j - 1, -1, -1):
            S = U[i, i] + U[j, j]
            if _np.isclose(S, 0):  # then use TU = UT
                if _np.isclose(T[i, i] - T[j, j], 0):  # then just set to zero
                    U[i, j] = 0.0  # TODO: check correctness of this case
                else:
                    U[i, j] = T[i, j] * (U[i, i] - U[j, j]) / (T[i, i] - T[j, j]) + \
                        sum([U[i, k] * T[k, j] - T[i, k] * U[k, j] for k in range(i + 1, j)]) \
                        / (T[i, i] - T[j, j])
            else:  # use U^2 = I
                U[i, j] = - sum([U[i, k] * U[k, j] for k in range(i + 1, j)]) / S
    return _np.dot(Z, _np.dot(U, _np.conjugate(Z.T)))

    #Quick & dirty - not always stable:
    #U,_,Vt = _np.linalg.svd(M)
    #return _np.dot(U,Vt) 
Example 15
def fun(self, x):
        dx, dy, dz = self._compute_coordinate_deltas(x)
        with np.errstate(divide='ignore'):
            dm1 = (dx**2 + dy**2 + dz**2) ** -0.5
        dm1[np.diag_indices_from(dm1)] = 0
        return 0.5 * np.sum(dm1) 
Example 16
def grad(self, x):
        dx, dy, dz = self._compute_coordinate_deltas(x)

        with np.errstate(divide='ignore'):
            dm3 = (dx**2 + dy**2 + dz**2) ** -1.5
        dm3[np.diag_indices_from(dm3)] = 0

        grad_x = -np.sum(dx * dm3, axis=1)
        grad_y = -np.sum(dy * dm3, axis=1)
        grad_z = -np.sum(dz * dm3, axis=1)

        return np.hstack((grad_x, grad_y, grad_z)) 
Example 17
Project: Mastering-Elasticsearch-7.0   Author: PacktPublishing   File: test_pairwise.py    License: MIT License 5 votes vote down vote up
def test_cosine_distances():
    # Check the pairwise Cosine distances computation
    rng = np.random.RandomState(1337)
    x = np.abs(rng.rand(910))
    XA = np.vstack([x, x])
    D = cosine_distances(XA)
    assert_array_almost_equal(D, [[0., 0.], [0., 0.]])
    # check that all elements are in [0, 2]
    assert np.all(D >= 0.)
    assert np.all(D <= 2.)
    # check that diagonal elements are equal to 0
    assert_array_almost_equal(D[np.diag_indices_from(D)], [0., 0.])

    XB = np.vstack([x, -x])
    D2 = cosine_distances(XB)
    # check that all elements are in [0, 2]
    assert np.all(D2 >= 0.)
    assert np.all(D2 <= 2.)
    # check that diagonal elements are equal to 0 and non diagonal to 2
    assert_array_almost_equal(D2, [[0., 2.], [2., 0.]])

    # check large random matrix
    X = np.abs(rng.rand(1000, 5000))
    D = cosine_distances(X)
    # check that diagonal elements are equal to 0
    assert_array_almost_equal(D[np.diag_indices_from(D)], [0.] * D.shape[0])
    assert np.all(D >= 0.)
    assert np.all(D <= 2.) 
Example 18
Project: Mastering-Elasticsearch-7.0   Author: PacktPublishing   File: pairwise.py    License: MIT License 5 votes vote down vote up
def cosine_distances(X, Y=None):
    """Compute cosine distance between samples in X and Y.

    Cosine distance is defined as 1.0 minus the cosine similarity.

    Read more in the :ref:`User Guide <metrics>`.

    Parameters
    ----------
    X : array_like, sparse matrix
        with shape (n_samples_X, n_features).

    Y : array_like, sparse matrix (optional)
        with shape (n_samples_Y, n_features).

    Returns
    -------
    distance matrix : array
        An array with shape (n_samples_X, n_samples_Y).

    See also
    --------
    sklearn.metrics.pairwise.cosine_similarity
    scipy.spatial.distance.cosine : dense matrices only
    """
    # 1.0 - cosine_similarity(X, Y) without copy
    S = cosine_similarity(X, Y)
    S *= -1
    S += 1
    np.clip(S, 0, 2, out=S)
    if X is Y or Y is None:
        # Ensure that distances between vectors and themselves are set to 0.0.
        # This may not be the case due to floating point rounding errors.
        S[np.diag_indices_from(S)] = 0.0
    return S


# Paired distances 
Example 19
Project: pyGeoStatistics   Author: whimian   File: krige2d.py    License: MIT License 5 votes vote down vote up
def block_covariance(self):
        "the block covariance"
        if self._block_covariance is None:
            self._block_covariance = 0
            if self.ndb <= 1:  # point kriging
                self._block_covariance = self.unbias
            else:  # block kriging
                cov = list()
                for x1, y1 in zip(self.xdb, self.ydb):
                    for x2, y2 in zip(self.xdb, self.ydb):
                        cov.append(self._cova2(x1, y1, x2, y2))
                cov = np.array(cov).reshape((self.ndb, self.ndb))
                cov[np.diag_indices_from(cov)] -= self.c0
                self._block_covariance = np.mean(cov)
        return self._block_covariance 
Example 20
Project: pyGeoStatistics   Author: whimian   File: cokrige.py    License: MIT License 5 votes vote down vote up
def block_covariance(self):
        "return average covariance within block"
        if self._block_covariance is None:
            if self.ndb <= 1:  # point kriging
                self._block_covariance = self.unbias
            else:
                cov = list()
                for x1, y1, z1 in zip(self.xdb, self.ydb, self.zdb):
                    for x2, y2, z2 in zip(self.xdb, self.ydb, self.zdb):
                        cov.append(self._cova3((x1, y1, z1), (x2, y2, z2)))
                cov = np.array(cov).reshape((self.ndb, self.ndb))
                cov[np.diag_indices_from(cov)] -= self.c0
                self._block_covariance = np.mean(cov)
        return self._block_covariance 
Example 21
Project: kite   Author: pyrocko   File: test_util.py    License: GNU General Public License v3.0 5 votes vote down vote up
def test_trim_matrix():
    arr = num.full((100, 100), num.nan)

    arr[-1, -1] = 1.
    assert util.trimMatrix(arr).shape == (1, 1)

    arr[-2, -2] = 1.
    assert util.trimMatrix(arr).shape == (2, 2)

    arr[num.diag_indices_from(arr)] = 1.
    assert util.trimMatrix(arr).shape == arr.shape 
Example 22
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_minimize_constrained.py    License: MIT License 5 votes vote down vote up
def fun(self, x):
        dx, dy, dz = self._compute_coordinate_deltas(x)
        with np.errstate(divide='ignore'):
            dm1 = (dx**2 + dy**2 + dz**2) ** -0.5
        dm1[np.diag_indices_from(dm1)] = 0
        return 0.5 * np.sum(dm1) 
Example 23
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_minimize_constrained.py    License: MIT License 5 votes vote down vote up
def grad(self, x):
        dx, dy, dz = self._compute_coordinate_deltas(x)

        with np.errstate(divide='ignore'):
            dm3 = (dx**2 + dy**2 + dz**2) ** -1.5
        dm3[np.diag_indices_from(dm3)] = 0

        grad_x = -np.sum(dx * dm3, axis=1)
        grad_y = -np.sum(dy * dm3, axis=1)
        grad_z = -np.sum(dz * dm3, axis=1)

        return np.hstack((grad_x, grad_y, grad_z)) 
Example 24
Project: nni   Author: microsoft   File: bayesian.py    License: MIT License 5 votes vote down vote up
def incremental_fit(self, train_x, train_y):
        """ Incrementally fit the regressor. """
        if not self._first_fitted:
            raise ValueError(
                "The first_fit function needs to be called first.")

        train_x, train_y = np.array(train_x), np.array(train_y)

        # Incrementally compute K
        up_right_k = edit_distance_matrix(self._x, train_x)
        down_left_k = np.transpose(up_right_k)
        down_right_k = edit_distance_matrix(train_x)
        up_k = np.concatenate((self._distance_matrix, up_right_k), axis=1)
        down_k = np.concatenate((down_left_k, down_right_k), axis=1)
        temp_distance_matrix = np.concatenate((up_k, down_k), axis=0)
        k_matrix = bourgain_embedding_matrix(temp_distance_matrix)
        diagonal = np.diag_indices_from(k_matrix)
        diagonal = (diagonal[0][-len(train_x):], diagonal[1][-len(train_x):])
        k_matrix[diagonal] += self.alpha

        try:
            self._l_matrix = cholesky(k_matrix, lower=True)  # Line 2
        except LinAlgError:
            return self

        self._x = np.concatenate((self._x, train_x), axis=0)
        self._y = np.concatenate((self._y, train_y), axis=0)
        self._distance_matrix = temp_distance_matrix

        self._alpha_vector = cho_solve(
            (self._l_matrix, True), self._y)  # Line 3

        return self 
Example 25
Project: qml   Author: qmlcode   File: math.py    License: MIT License 5 votes vote down vote up
def cho_solve(A, y):
    """ Solves the equation

            :math:`A x = y`

        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module. Preserves the input matrix A.

        :param A: Matrix (symmetric and positive definite, left-hand side).
        :type A: numpy array
        :param y: Vector (right-hand side of the equation).
        :type y: numpy array

        :return: The solution vector.
        :rtype: numpy array
        """

    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
        raise ValueError('expected square matrix')

    if len(y.shape) != 1 or y.shape[0] != A.shape[1]:
        raise ValueError('expected matrix and vector of same stride size')

    n = A.shape[0]

    # Backup diagonal before Cholesky-decomposition
    A_diag = A[np.diag_indices_from(A)]

    x = np.zeros(n)
    fcho_solve(A, y, x)

    # Reset diagonal after Cholesky-decomposition
    A[np.diag_indices_from(A)] = A_diag

    # Copy lower triangle to upper
    i_lower = np.tril_indices_from(A)
    A.T[i_lower] = A[i_lower]

    return x 
Example 26
Project: qml   Author: qmlcode   File: math.py    License: MIT License 5 votes vote down vote up
def bkf_solve(A, y):
    """ Solves the equation

            :math:`A x = y`

        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module. Preserves the input matrix A.

        :param A: Matrix (symmetric and positive definite, left-hand side).
        :type A: numpy array
        :param y: Vector (right-hand side of the equation).
        :type y: numpy array

        :return: The solution vector.
        :rtype: numpy array
        """

    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
        raise ValueError('expected square matrix')

    if len(y.shape) != 1 or y.shape[0] != A.shape[1]:
        raise ValueError('expected matrix and vector of same stride size')

    n = A.shape[0]

    # Backup diagonal before Cholesky-decomposition
    A_diag = A[np.diag_indices_from(A)]

    x = np.zeros(n)
    fbkf_solve(A, y, x)

    # Reset diagonal after Cholesky-decomposition
    A[np.diag_indices_from(A)] = A_diag

    # Copy lower triangle to upper
    i_lower = np.tril_indices_from(A)
    A.T[i_lower] = A[i_lower]

    return x 
Example 27
Project: Splunking-Crime   Author: nccgroup   File: pairwise.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def cosine_distances(X, Y=None):
    """Compute cosine distance between samples in X and Y.

    Cosine distance is defined as 1.0 minus the cosine similarity.

    Read more in the :ref:`User Guide <metrics>`.

    Parameters
    ----------
    X : array_like, sparse matrix
        with shape (n_samples_X, n_features).

    Y : array_like, sparse matrix (optional)
        with shape (n_samples_Y, n_features).

    Returns
    -------
    distance matrix : array
        An array with shape (n_samples_X, n_samples_Y).

    See also
    --------
    sklearn.metrics.pairwise.cosine_similarity
    scipy.spatial.distance.cosine (dense matrices only)
    """
    # 1.0 - cosine_similarity(X, Y) without copy
    S = cosine_similarity(X, Y)
    S *= -1
    S += 1
    np.clip(S, 0, 2, out=S)
    if X is Y or Y is None:
        # Ensure that distances between vectors and themselves are set to 0.0.
        # This may not be the case due to floating point rounding errors.
        S[np.diag_indices_from(S)] = 0.0
    return S


# Paired distances 
Example 28
Project: tape   Author: songlab-cal   File: datasets.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def __getitem__(self, index):
        item = self.data[index]

        msa = item['msa']
        dist = item['dist6d']
        omega = item['omega6d']
        theta = item['theta6d']
        phi = item['phi6d']

        if self._split == 'train':
            msa = self._subsample_msa(msa)
        elif self._split == 'valid':
            msa = msa[:20000]  # runs out of memory if msa is way too big
        msa, dist, omega, theta, phi = self._slice_long_sequences(
            msa, dist, omega, theta, phi)

        mask = dist == 0

        dist_bins = np.digitize(dist, self._dist_bins)
        omega_bins = np.digitize(omega, self._dihedral_bins) + 1
        theta_bins = np.digitize(theta, self._dihedral_bins) + 1
        phi_bins = np.digitize(phi, self._planar_bins) + 1

        dist_bins[mask] = 0
        omega_bins[mask] = 0
        theta_bins[mask] = 0
        phi_bins[mask] = 0

        dist_bins[np.diag_indices_from(dist_bins)] = -1

        # input_mask = np.ones_like(msa[0])

        return msa, dist_bins, omega_bins, theta_bins, phi_bins 
Example 29
Project: snfpy   Author: rmarkello   File: compute.py    License: GNU Lesser General Public License v3.0 5 votes vote down vote up
def get_n_clusters(arr, n_clusters=range(2, 6)):
    """
    Finds optimal number of clusters in `arr` via eigengap method

    Parameters
    ----------
    arr : (N, N) array_like
        Input array (e.g., the output of :py:func`snf.compute.snf`)
    n_clusters : array_like
        Numbers of clusters to choose between

    Returns
    -------
    opt_cluster : int
        Optimal number of clusters
    second_opt_cluster : int
        Second best number of clusters
    """

    # confirm inputs are appropriate
    n_clusters = check_array(n_clusters, ensure_2d=False)
    n_clusters = n_clusters[n_clusters > 1]

    # don't overwrite provided array!
    graph = arr.copy()

    graph = (graph + graph.T) / 2
    graph[np.diag_indices_from(graph)] = 0
    degree = graph.sum(axis=1)
    degree[np.isclose(degree, 0)] += np.spacing(1)
    di = np.diag(1 / np.sqrt(degree))
    laplacian = di @ (np.diag(degree) - graph) @ di

    # perform eigendecomposition and find eigengap
    eigs = np.sort(np.linalg.eig(laplacian)[0])
    eigengap = np.abs(np.diff(eigs))
    eigengap = eigengap * (1 - eigs[:-1]) / (1 - eigs[1:])
    n = eigengap[n_clusters - 1].argsort()[::-1]

    return n_clusters[n[:2]] 
Example 30
Project: domain_adversarial_neural_network   Author: GRAAL-Research   File: mSDA.py    License: BSD 2-Clause "Simplified" License 5 votes vote down vote up
def mda_fit(X, noise=0.5, eta=1e-5):
    """
    inputs: 
        X : d x n input (Transpose of the usual data-matrix)
        noise: corruption level
        eta: regularization 
    
    outputs:
        hx: d x n hidden representation
        W: d x (d+1) mapping
    """
    d, n = np.shape(X)
    
    # adding bias
    Xb = np.vstack((X, np.ones(n)))
    
    # scatter matrix S
    S = np.dot(Xb, Xb.T)
    
    # corruption vector
    q = np.ones((d+1, 1)) * (1.-noise)
    q[-1] = 1
    
    # Q: (d+1)x(d+1)
    Q = S*np.dot(q,q.T)
    Q[np.diag_indices_from(Q)] = q.T[0] * np.diag(S)

    #P: dx(d+1)
    P = S[0:-1,:] * q.T 
    
    # final W = P*Q^-1, dx(d+1)
    reg = eta * np.eye(d+1)
    reg[-1,-1] = 0
    W = np.linalg.solve(Q.T + reg, P.T).T

    hx = np.tanh(np.dot(W, Xb))
    return hx, W