Python numpy.diag_indices_from() 使用实例

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 f_electric(P):
    n = len(P)
    a = 1e-9*np.array([p.a for p in P])
    apa = a[:, None] + a[None, :]


    x = 1e-9*np.array([p.x.flatten() for p in P])

    R = x[:, None, :] - x[None, :, :]
    r = np.sqrt(np.sum(R**2, 2) + 1e-100)
    R0 = R / (r**3)[:, :, None]

    q = np.array([float(p.charge) for p in P])
    const = qq**2 / (4.*np.pi*eperm*rpermw)

    QQ = q[:, None] * q[None, :]
    F = const * QQ[:, :, None] * R0
    #F[np.diag_indices_from(r)] = 0.
    tooclose = r <= apa
    R0i = R / (np.maximum(a[:, None], a[None, :])**3)[:, :, None]
    F[tooclose] = (const * QQ[:, :, None] * R0i)[tooclose]

    f = np.sum(F, 1).T.reshape(3*n, 1)

    return f 

Example 2

def kth_diag_indices(k, a):
	""" Get diagonal indices of 2D array 'a' offset by 'k'

	Parameters
	----------
	k : int
		Diagonal offset

	a : numpy.ndarray
		Input numpy 2D array

	Returns
	-------
	indices : tuple of two numpy.ndarray
		It contain indences of elements that are at the diagonal offset by 'k'.

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

Example 3

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 4

def _solve_hessian(G, Y, thY, precon, lambda_min):
    N, T = Y.shape
    # Compute the derivative of the score
    psidY = ne.evaluate('(- thY ** 2 + 1.) / 2.')  # noqa
    # Build the diagonal of the Hessian, a.
    Y_squared = Y ** 2
    if precon == 2:
        a = np.inner(psidY, Y_squared) / float(T)
    elif precon == 1:
        sigma2 = np.mean(Y_squared, axis=1)
        psidY_mean = np.mean(psidY, axis=1)
        a = psidY_mean[:, None] * sigma2[None, :]
        diagonal_term = np.mean(Y_squared * psidY) + 1.
        a[np.diag_indices_from(a)] = diagonal_term
    else:
        raise ValueError('precon should be 1 or 2')
    # Compute the eigenvalues of the Hessian
    eigenvalues = 0.5 * (a + a.T - np.sqrt((a - a.T) ** 2 + 4.))
    # Regularize
    problematic_locs = eigenvalues < lambda_min
    np.fill_diagonal(problematic_locs, False)
    i_pb, j_pb = np.where(problematic_locs)
    a[i_pb, j_pb] += lambda_min - eigenvalues[i_pb, j_pb]
    # Invert the transform
    return (G * a.T - G.T) / (a * a.T - 1.) 

Example 5

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 izip(self.xdb, self.ydb, self.zdb):
                    for x2, y2, z2 in izip(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 6

def run(self, verbose=True):
        """
        Conducts particle swarm optimization

        :param verbose: indicates whether or not to print progress regularly
        :return: best member of swarm and objective function value of best member of swarm
        """
        self._clear()
        for i in range(self.max_steps):
            self.cur_steps += 1

            if verbose and ((i + 1) % 100 == 0):
                print(self)

            u1 = zeros((self.swarm_size, self.swarm_size))
            u1[diag_indices_from(u1)] = [random() for x in range(self.swarm_size)]
            u2 = zeros((self.swarm_size, self.swarm_size))
            u2[diag_indices_from(u2)] = [random() for x in range(self.swarm_size)]

            vel_new = (self.c1 * self.vel) + \
                      (self.c2 * dot(u1, (self.best - self.pos))) + \
                      (self.c3 * dot(u2, (self.global_best - self.pos)))

            pos_new = self.pos + vel_new

            self._best(self.pos, pos_new)
            self.pos = pos_new
            self.scores = self._score(self.pos)
            self._global_best()

            if self._objective(self.global_best[0]) < (self.min_objective or 0):
                print("TERMINATING - REACHED MINIMUM OBJECTIVE")
                return self.global_best[0], self._objective(self.global_best[0])
        print("TERMINATING - REACHED MAXIMUM STEPS")
        return self.global_best[0], self._objective(self.global_best[0]) 

Example 7

def build_base_operator(self, t):
        """
        :param t: Not used as mu and sigma are constant
        :return:
        """
        # Update drift and volatility
        self.build_moment_vectors(t)

        base_operator = np.zeros((self.d, self.d))

        nabla = linalg.block_diag(*[self.build_gradient_matrix(x) for x in range(1, self.d - 1)])

        moments = np.zeros(2 * (self.d - 2))
        for i in range(0, self.d - 2):
            moments[2 * i] = self.drift[i + 1]
            moments[2 * i + 1] = self.volatility[i + 1]

        generator_elements = linalg.solve(nabla, moments)

        r_idx, c_idx = np.diag_indices_from(base_operator)
        base_operator[r_idx[1:-1], c_idx[:-2]] = generator_elements[::2]
        base_operator[r_idx[1:-1], c_idx[2:]] = generator_elements[1::2]
        np.fill_diagonal(base_operator, -np.sum(base_operator, axis=1))

        # -- Boundary Condition: Volatility Matching --
        nabla_0 = self.grid[1] - self.grid[0]
        base_operator[0, 0] = - 1. * self.volatility[0] / (nabla_0 * nabla_0)
        base_operator[0, 1] = - base_operator[0, 0]

        nabla_d = self.grid[self.d - 1] - self.grid[self.d - 2]
        base_operator[self.d - 1, self.d - 1] = - 1. * self.volatility[self.d - 1] / (nabla_d * nabla_d)
        base_operator[self.d - 1, self.d - 2] = - base_operator[self.d - 1, self.d - 1]
        # ----------------------------------------------

        self.sanity_check_base_operator(base_operator)

        return base_operator 

Example 8

def solve(self):

		DI = np.diag_indices_from(self.SII)

		D = self.SII[DI]
		#for i in xrange(1, self.SII.shape[0]): self.SII[i:,i-1] = self.SII[i-1,i:]

		self.SII[DI] += np.mean(D)
		_, self.W, _ = sposv(self.SII, self.SIO, overwrite_a=True, overwrite_b=False) #, lower=1)

		#self.SII[DI] = D 

Example 9

def divide_diagonal_by_2(CHI0, div_fact = 2.):
    CHI = CHI0.copy();
    CHI[np.diag_indices_from(CHI)] /= div_fact
    return CHI 

Example 10

def _add_to_diag(A, z):
    B = A.copy()
    B[np.diag_indices_from(B)] + z
    return B 

Example 11

def from_quat(quat):
    """Create a direction cosine matrix from a quaternion.

    First 3 elements of the quaternion form its vector part.

    Parameters
    ----------
    quat : array_like, shape (4,) or (n, 4)
        Quaternions.

    Returns
    -------
    dcm : ndarray, shape (3, 3) or (n, 3, 3)
        Direction cosine matrices.
    """
    q = np.asarray(quat)
    if q.ndim == 1:
        rho = q[:3]
        q4 = q[3]
        rho_skew = _skew_matrix_single(rho)
        dcm = 2 * (np.outer(rho, rho) + q4 * rho_skew)
        dcm[np.diag_indices_from(dcm)] += q4**2 - np.dot(rho, rho)
    else:
        rho = q[:, :3]
        q4 = q[:, 3]
        rho_skew = _skew_matrix_array(rho)
        dcm = 2 * (rho[:, None, :] * rho[:, :, None] +
                   q4[:, None, None] * rho_skew)
        diag = q4**2 - np.sum(rho**2, axis=1)
        dcm[:, np.arange(3), np.arange(3)] += diag[:, None]

    return dcm 

Example 12

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 13

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 14

def tail_correlations(corrmat, tails, mask=None):
    """Compute the mean within and between tail correlations."""
    assert corrmat.shape == (len(tails), len(tails))
    if mask is not None:
        assert corrmat.shape == mask.shape

    def is_within(pair):
        a, b = pair
        return a == b

    # Identify cells in the matrix that represent within or between pairs
    pairs = product(tails, tails)
    wmat = np.reshape(map(is_within, pairs), corrmat.shape)
    bmat = ~wmat

    # Possibly exclude cells with the mask
    if mask is not None:
        wmat &= mask
        bmat &= mask

    # Remove the diagonal from the within matrix
    wmat[np.diag_indices_from(wmat)] = False

    # Average the correlations for within and between pairs
    corr_w = corrmat[wmat].mean()
    corr_b = corrmat[bmat].mean()

    return corr_w, corr_b 

Example 15

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 izip(self.xdb, self.ydb):
                    for x2, y2 in izip(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 16

def kth_diag_indices(array, diag_k):
    """Return a tuple of indices for retrieving the k'th diagonal
    of matrix a.

    :param array: Input matrix.
    :type array: :class:`numpy array <numpy.ndarray>`
    :param int diag_k: Diagonal to index. 0 is the centre, 1 is the first \
            diagonal below the centre, -1 is the first diagonal above \
            the centre.

    >>> my_matrix = np.array([[ 0, -1, -2, -3],
    ...                       [ 1,  0, -1, -2],
    ...                       [ 2,  1,  0, -1],
    ...                       [ 3,  2,  1,  0]])
    >>> matrix.kth_diag_indices(my_matrix, 1)
    (array([1, 2, 3]), array([0, 1, 2]))
    >>> my_matrix[
    ...    matrix.kth_diag_indices(my_matrix, 1)
    ...  ]
    array([1, 1, 1])
    >>> my_matrix[
    ...    matrix.kth_diag_indices(my_matrix, -2)
    ...  ]
    array([-2, -2])
    """

    rows, cols = np.diag_indices_from(array)
    if diag_k < 0:
        return rows[:diag_k], cols[-diag_k:]

    if diag_k > 0:
        return rows[diag_k:], cols[:-diag_k]

    return rows, cols 

Example 17

def evaluate(self, x):
		try:
			return multivariate_normal.logpdf(x, self.mean, self.cov)
		except:
			self.cov[np.diag_indices_from(self.cov)] += 1e-4
			return multivariate_normal.logpdf(x, self.mean, self.cov) 

Example 18

def _calculate_reduced_likelihood_params(self, thetas=None):
        """
        Calculate quantity with same maximum location as the log-likelihood for a given theta.

        Parameters
        ----------
        thetas : ndarray, optional
            Given input correlation coefficients. If none given, uses self.thetas
            from training.
        """
        if thetas is None:
            thetas = self.thetas

        X, Y = self.X, self.Y
        params = {}

        # Correlation Matrix
        distances = np.zeros((self.n_samples, self.n_dims, self.n_samples))
        for i in range(self.n_samples):
            distances[i, :, i + 1:] = np.abs(X[i, ...] - X[i + 1:, ...]).T
            distances[i + 1:, :, i] = distances[i, :, i + 1:].T

        R = np.exp(-thetas.dot(np.square(distances)))
        R[np.diag_indices_from(R)] = 1. + self.nugget

        [U, S, Vh] = linalg.svd(R)

        # Penrose-Moore Pseudo-Inverse:
        # Given A = USV^* and Ax=b, the least-squares solution is
        # x = V S^-1 U^* b.
        # Tikhonov regularization is used to make the solution significantly
        # more robust.
        h = 1e-8 * S[0]
        inv_factors = S / (S ** 2. + h ** 2.)

        alpha = Vh.T.dot(np.einsum('j,kj,kl->jl', inv_factors, U, Y))
        logdet = -np.sum(np.log(inv_factors))
        sigma2 = np.dot(Y.T, alpha).sum(axis=0) / self.n_samples
        reduced_likelihood = -(np.log(np.sum(sigma2)) +
                               logdet / self.n_samples)

        params['alpha'] = alpha
        params['sigma2'] = sigma2 * np.square(self.Y_std)
        params['S_inv'] = inv_factors
        params['U'] = U
        params['Vh'] = Vh

        return reduced_likelihood, params 

Example 19

def test_krr_cmat():

    test_dir = os.path.dirname(os.path.realpath(__file__))

    # Parse file containing PBE0/def2-TZVP heats of formation and xyz filenames
    data = get_energies(test_dir + "/data/hof_qm7.txt")

    # Generate a list of qml.Compound() objects
    mols = []

    for xyz_file in sorted(data.keys())[:1000]:

        # Initialize the qml.Compound() objects
        mol = qml.Compound(xyz=test_dir + "/qm7/" + xyz_file)

        # Associate a property (heat of formation) with the object
        mol.properties = data[xyz_file]

        # This is a Molecular Coulomb matrix sorted by row norm
        mol.generate_coulomb_matrix(size=23, sorting="row-norm")

        mols.append(mol)

    # Shuffle molecules
    np.random.seed(666)
    np.random.shuffle(mols)

    # Make training and test sets
    n_test  = 300
    n_train = 700

    training = mols[:n_train]
    test  = mols[-n_test:]

    # List of representations
    X  = np.array([mol.representation for mol in training])
    Xs = np.array([mol.representation for mol in test])

    # List of properties
    Y = np.array([mol.properties for mol in training])
    Ys = np.array([mol.properties for mol in test])

    # Set hyper-parameters
    sigma = 10**(4.2)
    llambda = 10**(-10.0)

    # Generate training Kernel
    K = laplacian_kernel(X, X, sigma)

    # Solve alpha
    K[np.diag_indices_from(K)] += llambda
    alpha = cho_solve(K,Y)

    # Calculate prediction kernel
    Ks = laplacian_kernel(X, Xs, sigma)
    Yss = np.dot(Ks.transpose(), alpha)

    mae = np.mean(np.abs(Ys - Yss))
    print(mae) 
点赞

发表评论

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