Python numpy 模块，diag_indices_from() 实例源码

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

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

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

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.)

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

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])

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

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

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

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

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

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)

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))

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

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

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

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)

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

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)