Python numpy.diagonal函数代码示例

def main():
    mu = np.array([[5],
                   [6],
                   [7],
                   [8],
                   ])
    S = np.array([[3.01602775, 1.02746769, -3.60224613, -2.08792829],
                  [1.02746769, 5.65146472, -3.98616664, 0.48723704],
                  [-3.60224613, -3.98616664, 13.04508284, -1.59255406],
                  [-2.08792829, 0.48723704, -1.59255406, 8.28742469],
                  ])
    mu_est, S_est = estimate(mu, S)

    print('mean estimation:', mu_est)
    print('S estimation:')
    print(S_est)

    print()

    print('mean % difference:', (np.diagonal(mu / mu_est) - 1) * 100)
    print('S % difference:')
    print(((S / S_est) - 1) * 100)

    print()

    mu_est, S_est = estimate_more(mu, S)
    print('mean estimation:', mu_est)
    print('S estimation:')
    print(S_est)

    print()

    print('mean % difference:', (np.diagonal(mu / mu_est) - 1) * 100)
    print('S % difference:')
    print(((S / S_est) - 1) * 100)

def checkSquare(board, playerChar):
    otherPlayerChar = notPlayer(playerChar)
    lookupBoard = numpy.array(board)
    foundLocation = (INVALID, INVALID)
    
    # Rows and cols
    for i in range(0, BOARD_SIZE):
        #print("        Processing square index " + str(i))
        rowSet = lookupBoard[i, :].tolist()
        foundIndex = checkSet(rowSet, otherPlayerChar)
        if (foundIndex != INVALID):
            foundLocation = (i, foundIndex)
        
        colSet = lookupBoard[:, i].tolist()
        foundIndex = checkSet(colSet, otherPlayerChar)
        if (foundIndex >= 0):
            foundLocation = (foundIndex, i)

    # Diagonals
    foundIndex = checkSet(numpy.diagonal(lookupBoard).tolist(), otherPlayerChar) 
    if (foundIndex != INVALID):
        foundLocation = (foundIndex, foundIndex)
        
    foundIndex = checkSet(numpy.diagonal(lookupBoard[::-1]).tolist(), otherPlayerChar) 
    if (foundIndex != INVALID):
        revFoundIndex = BOARD_SIZE - foundIndex - 1
        foundLocation = (revFoundIndex, foundIndex)
                        
    #print(foundLocation)
    return foundLocation

def kl_Divergence(mean1, cov1, mean2, cov2):
	N = mean1.size
	#check length of the mean vectors
	if N != mean2.size:
		raise Exception("Mean sizes do not match!")
	#check that cov matrices have the same length as the mean
	if cov1.shape[0] != N or cov1.shape[1] != N:
		raise Exception("cov1 sizes do not equal mean length!")
	if cov2.shape[0] != N or cov2.shape[1] != N:
		raise Exception("cov2 sizes do not equal mean length!")

	#return Cholesky decompositions for covariance matrices
	chol1 = np.linalg.cholesky(cov1)
	chol2 = np.linalg.cholesky(cov2)
	#begin distance calculation
	ld1 = 2 * np.sum(np.log10(np.diagonal(chol1)), axis=0)
	ld2 = 2 * np.sum(np.log10(np.diagonal(chol2)), axis=0)

	#calculate det
	ldet = ld2 - ld1
	#inverse from Cholesky decomposition
	S1i = np.dot((np.linalg.inv(np.transpose(chol1))), np.linalg.inv(chol1))
	tr = np.sum(np.diagonal(np.dot(S1i, cov2)), axis=0)
	m2mm1 = np.subtract(mean2, mean1)
	
	#asNumeric equivalent in python...
	qf = np.dot(np.transpose(m2mm1), np.dot(S1i, m2mm1))
	r = 0.5 * (ldet + tr + qf - N)
	return r

def plot_classifier(cal,trues,falses):
    x = np.linspace(-1,1,200)
    y = np.linspace(-1,1,200)
    (X,Y) = np.meshgrid(x,y)
    XY = np.dstack([X,Y])
    Z = lr.evaluate_poly(cal,lr.dim2_deg4[:,None,None,:],XY)

    f = plt.figure()
    ax = f.add_subplot(111)
    ax.grid(True)
    ax.axis('equal')

    ax.contour(X,Y,Z,levels=np.linspace(np.min(Z),np.max(Z),20),colors='k')

    # adjust for the way the random samples were distributed
    adjustment = 1

    if len(falses)>0:
        falses = np.diagonal(falses,axis1=1,axis2=2)[:,:2] - [adjustment]*2
        ax.scatter(*falses.T,color='red')

    if len(trues)>0:
        trues = np.diagonal(trues,axis1=1,axis2=2)[:,:2] - [adjustment]*2
        ax.scatter(*trues.T,color='green')

    return f

    def board_status(b):
        '''
            Checks for a winner of the game.
        '''
        board = np.array(b)
        for r in range(3):
            row = board[r]
            column = board[:,r]
            if np.all(row == user.X):
                return user.X
            if np.all(row == user.O):
                return user.O
            if np.all(column == user.X):
                return user.X
            if np.all(column == user.O):
                return user.O

        diagonal = np.diagonal(board)

        if np.all(diagonal == user.O):
            return user.O

        if np.all(diagonal == user.X):
            return user.X

        board_flip =  np.fliplr(board)
        reverse_diagonal = np.diagonal(board_flip)

        if np.all(reverse_diagonal == user.O):
            return user.O

        if np.all(reverse_diagonal == user.X):
            return user.X

        return user.available

def test():
    # feed = DataFeed("/mnt/e/data/algae_dataset", "/mnt/e/data/algae_patches", False)
    # feed = DataFeed("/mnt/e/data/algae_dataset_cells_only", "/mnt/e/data/algae_patches_cells_only", False, False)
    feed = EqualizedDataFeed("/mnt/e/data/algae_dataset_equal_batches", False)
    sess = tf.InteractiveSession()
    batch_size = 200
    tile_size = (128,128)
    lr = 1e-3
    eps = 1.0
    a = 0.001
    net = AlexNet(batch_size, tile_size, sess, lr, eps, a)
    saver = tf.train.Saver(net.train_vars)

    saver.restore(sess, "AlexNet_4ft_1mp_2fc_1softmax_random_symmetry_201611141618_best.ckpt")

    conf_mat = np.zeros((4,4))

    for i in range(5):
        X,Y = feed.next_batch(batch_size)
        pred = net.forward(X).argmax(axis=1)
        
        targets = Y.argmax(axis=1)
        
        for i in range(4):
            for j in range(4):
                conf_mat[i,j] += ((pred==i)*(targets==j)).sum()

    print conf_mat
    print np.diagonal(conf_mat).sum()*1./(conf_mat.sum())

  def compute_distances_no_loops(self, X):
    """
    Compute the distance between each test point in X and each training point
    in self.X_train using no explicit loops.

    Input / Output: Same as compute_distances_two_loops
    """
    num_test = X.shape[0]
    num_train = self.X_train.shape[0]
    dists = np.zeros((num_test, num_train)) 
    #########################################################################
    # TODO:                                                                 #
    # Compute the l2 distance between all test points and all training      #
    # points without using any explicit loops, and store the result in      #
    # dists.                                                                #
    #                                                                       #
    # You should implement this function using only basic array operations; #
    # in particular you should not use functions from scipy.                #
    #                                                                       #
    # HINT: Try to formulate the l2 distance using matrix multiplication    #
    #  and two broadcast sums.                                         #
    #########################################################################
    dists = -2*np.dot(X, np.transpose(self.X_train))
    assert(dists.shape == (num_test, num_train))
    train_diags = np.diagonal(np.dot(self.X_train, np.transpose(self.X_train)))
    assert(train_diags.shape == (num_train,))
    test_diags = np.diagonal(np.dot(X, np.transpose(X)))
    dists = dists + train_diags # broadcast sum #1
    dists = np.transpose(np.transpose(dists) + test_diags)  # broadcast sum #2
    #########################################################################
    #                         END OF YOUR CODE                              #
    #########################################################################
    return np.sqrt(dists)

def main():

    sample='q'
    sm_bin='10.0_10.5'
    catalogue = 'sm_9.5_s0.2_sfr_c-0.75_250'

    #load in fiducial mock
    filepath = './'
    filename = 'sm_9.5_s0.2_sfr_c-0.8_Chinchilla_250_wp_fiducial_'+sample+'_'+sm_bin+'_cov.npy'
    cov = np.matrix(np.load(filepath+filename))
    diag = np.diagonal(cov)
    filepath = cu.get_output_path() + 'analysis/central_quenching/observables/'
    filename = 'sm_9.5_s0.2_sfr_c-0.8_Chinchilla_250_wp_fiducial_'+sample+'_'+sm_bin+'.dat'
    data = ascii.read(filepath+filename)
    rbins = np.array(data['r'])
    mu = np.array(data['wp'])
    
    #load in comparison mock
    
    
    
    
    plt.figure()
    plt.errorbar(rbins, mu, yerr=np.sqrt(np.diagonal(cov)), color='black')
    plt.plot(rbins, wp,  color='red')
    plt.xscale('log')
    plt.yscale('log')
    plt.show()
    
    inv_cov = cov.I
    Y = np.matrix((wp-mu))
    
    X = Y*inv_cov*Y.T
    
    print(X)

    def getWorldSIP(self, position):
        u = position[0] - self.pixReference[0]
        v = position[1] - self.pixReference[1]
        fuvArray = numpy.zeros((4, 4))
        fuvArray[0][0] = 1
        fuvArray[0][1] = v
        fuvArray[0][2] = v * v
        fuvArray[0][3] = v * v * v
        fuvArray[1][0] = u
        fuvArray[1][1] = u * v
        fuvArray[1][2] = u * v * v
        fuvArray[2][0] = u * u
        fuvArray[2][1] = u * u * v
        fuvArray[3][0] = u * u * u

        uprime = numpy.sum(numpy.diagonal(numpy.dot(self.SIP_A, fuvArray)))
        vprime = numpy.sum(numpy.diagonal(numpy.dot(self.SIP_B, fuvArray)))

        # print "u': " + str(uprime)
        # print "v': " + str(vprime)
        u = u + uprime
        v = v + vprime
        CD = numpy.array(self.linearTransform)
        pixel = numpy.array((u, v))
        world = numpy.dot(CD, pixel)
        world = (world[0] + self.raDeg, world[1] + self.dec)

        return (world[0], world[1])

  def compute_distances_no_loops(self, X):
    """
    Compute the distance between each test point in X and each training point
    in self.X_train using no explicit loops.

    Input / Output: Same as compute_distances_two_loops
    """
    num_test = X.shape[0]
    num_train = self.X_train.shape[0]
    dists = np.zeros((num_test, num_train)) 
    X1=X
    X2=self.X_train
    #########################################################################
    # TODO:                                                                 #
    # Compute the l2 distance between all test points and all training      #
    # points without using any explicit loops, and store the result in      #
    # dists.                                                                #
    #                                                                       #
    # You should implement this function using only basic array operations; #
    # in particular you should not use functions from scipy.                #
    #                                                                       #
    # HINT: Try to formulate the l2 distance using matrix multiplication    #
    #       and two broadcast sums.                                         #
    #########################################################################
    X1X1Tdiag = np.array(np.diagonal(np.dot(X1,X1.T)))[np.newaxis]
    X1X1T = np.tile(X1X1Tdiag.T,(1,num_train))
    X2X2Tdiag = np.array(np.diagonal(np.dot(X2,X2.T)))[np.newaxis]
    X2X2T = np.tile(X2X2Tdiag,(num_test,1))
    X1X2T = np.dot(X1,X2.T)
    dists = np.sqrt(X1X1T - 2*X1X2T + X2X2T)

    #########################################################################
    #                         END OF YOUR CODE                              #
    #########################################################################
    return dists

    def mixing_constants(self, T, x):
        """ Compute Am, Bm, dAmdT and d2AmdT2,
        the "mixing constants" of the model and
        associated derivatives."""

        Tcmat = self.Tcmat
        Pcmat = self.Pcmat
        Vcmat = self.Vcmat
        Omat = self.Omat
        C = 0.37464 + 1.54226*Omat - 0.26992*Omat**2.0
        B = 0.07796*RGAS*(np.diagonal(Tcmat)/
                np.diagonal(Pcmat))
        A = (0.457236*((RGAS*Tcmat)**2.0)/Pcmat)*\
                ((1.0 +  C*(1.0 -
                    np.sqrt(np.tensordot(T, 1./Tcmat, 0))))**2.0)

        Am = np.tensordot((A*x[:,:,:,:,None]*x[:,:,:,None,:]),
                np.ones(Tcmat.shape), 2)
        Bm = np.tensordot(x, B, 1)

        G = C*np.sqrt(np.tensordot(T, 1./Tcmat, axes=0))/(
                1.0 + C*(1.0 -
                np.sqrt(np.tensordot(T, 1./Tcmat, axes=0))))
        
        dAmdT = (-1./T)*(np.tensordot(
            G*A*x[:,:,:,None,:]*x[:,:,:,:,None],
            np.ones(Tcmat.shape),
                    2))
        d2AmdT2 = 0.457236*(RGAS**2)/(2.0*T*np.sqrt(T))*\
            np.tensordot(
                (C*(1.+C)*Tcmat*np.sqrt(Tcmat)/Pcmat)*
                x[:,:,:,None,:]*x[:,:,:,:,None],
                np.ones(Tcmat.shape),
                                2)
        return Am, Bm, dAmdT, d2AmdT2

    def getPixel(self, position):
        pixel = (0, 0)
        u = position[0] - self.raDeg
        v = position[1] - self.dec
        # print "Incoming", u, v
        world = numpy.array((u, v))
        CD = numpy.array(self.linearTransform)
        invCD = numpy.linalg.inv(CD)
        pixel = numpy.dot(invCD, world)

        fuvArray = numpy.zeros((4, 4))
        fuvArray[0][0] = 1
        fuvArray[0][1] = v
        fuvArray[0][2] = v * v
        fuvArray[0][3] = v * v * v
        fuvArray[1][0] = u
        fuvArray[1][1] = u * v
        fuvArray[1][2] = u * v * v
        fuvArray[2][0] = u * u
        fuvArray[2][1] = u * u * v
        fuvArray[3][0] = u * u * u

        pixel[0] = pixel[0] + numpy.sum(numpy.diagonal(numpy.dot(self.SIP_AP, fuvArray)))
        pixel[1] = pixel[1] + numpy.sum(numpy.diagonal(numpy.dot(self.SIP_BP, fuvArray)))

        # print "Reference", self.pixReference, "offset", pixel
        pixel = pixel + self.pixReference
        return (pixel[0], pixel[1])

    def plot(self, eta1, u1, v1, eta2=None, u2=None, v2=None):
        self.fig.add_subplot(self.gs[0, 0])
        self.sp_eta.set_data(eta1)

        self.fig.add_subplot(self.gs[0, 1])
        self.sp_u.set_data(u1)

        self.fig.add_subplot(self.gs[0, 2])
        self.sp_v.set_data(v1)
            
        self.fig.add_subplot(self.gs[1, 0])
        self.sp_radial1.set_ydata(eta1.ravel());

        self.fig.add_subplot(self.gs[1, 1])
        self.sp_x_axis1.set_ydata(eta1[(self.ny+2)/2,:])
        self.sp_y_axis1.set_ydata(eta1[:,(self.nx+2)/2])

        self.fig.add_subplot(self.gs[1, 2])
        self.sp_x_diag1.set_ydata(np.diagonal(eta1, offset=-abs(self.nx-self.ny)/2))
        self.sp_y_diag1.set_ydata(np.diagonal(eta1.T, offset=abs(self.nx-self.ny)/2))
        
        if (eta2 is not None):
            self.fig.add_subplot(self.gs[2, 0])
            self.sp_radial2.set_ydata(eta2.ravel());

            self.fig.add_subplot(self.gs[2, 1])
            self.sp_x_axis2.set_ydata(eta2[(self.ny+2)/2,:])
            self.sp_y_axis2.set_ydata(eta2[:,(self.nx+2)/2])

            self.fig.add_subplot(self.gs[2, 2])
            self.sp_x_diag2.set_ydata(np.diagonal(eta2, offset=-abs(self.nx-self.ny)/2))
            self.sp_y_diag2.set_ydata(np.diagonal(eta2.T, offset=abs(self.nx-self.ny)/2))
        
        plt.draw()
        time.sleep(0.001)

    def __cal_shift_integrand(self, alpha=0, beta=0, gamma=0):
        """
        calculate shift current integrand and store it in 'shift_integrand'
        all parameters in this function are in hamiltonian gauge
        :param alpha, beta, gamma: 0: x, 1: y, 2: z
        """
        fermi_energy = self.fermi_energy
        nkpts = self.nkpts
        self.calculate('eigenvalue')
        self.calculate('A_h_ind', alpha)
        self.calculate('A_h_ind', beta)
        self.calculate('A_h_ind', gamma)
        self.calculate('A_h_ind_ind', beta, alpha)
        self.calculate('A_h_ind_ind', gamma, alpha)
        for i in range(nkpts):
            A_alpha = self.kpt_data['A_h_ind'][alpha][:, :, i]
            A_beta = self.kpt_data['A_h_ind'][beta][:, :, i]
            A_gamma = self.kpt_data['A_h_ind'][gamma][:, :, i]
            A_beta_alpha = self.kpt_data['A_h_ind_ind'][beta][alpha][:, :, i]
            A_gamma_alpha = self.kpt_data['A_h_ind_ind'][gamma][alpha][:, :, i]
            fermi = np.zeros(self.num_wann, dtype='float')
            fermi[self.kpt_data['eigenvalue'][:, i] > fermi_energy] = 0
            fermi[self.kpt_data['eigenvalue'][:, i] < fermi_energy] = 1
            fermi = fermi[:, None] - fermi[None, :]
            ki = np.diagonal(A_alpha)[:, None] - np.diagonal(A_alpha)[None, :]

            self.kpt_data['shift_integrand'][alpha][beta][gamma][:, :, i] = \
                np.imag(fermi *
                        (A_beta.T * (A_gamma_alpha - 1j * ki * A_gamma) +
                         A_gamma.T * (A_beta_alpha - 1j * ki * A_beta))
                        ) / 2

def gp_plot_prediction(
        predict_x,
        mean,
        variance=None
):
    """
    Plot a gp's prediction using pylab including error bars if variance specified

    Error bars are 2 * standard_deviation as in `Gaussian Processes for Machine Learning`__ by Rasmussen and Williams. 

    __ http://www.amazon.co.uk/Gaussian-Processes-Learning-Adaptive-Computation/dp/026218253X/
    """
    from pylab import plot, concatenate, fill
    if None != variance:
        # check variances are just about +ve - could signify a bug if not
        assert diagonal(variance).all() > -1e-10
        data = [
            (x, y, max(v, 0.0))
            for x, y, v
            in zip(predict_x, mean.flat, diagonal(variance))
        ]
    else:
        data = [
            (x, y)
            for x, y
            in zip(predict_x, mean)
        ]
    data.sort(key=lambda d: d[0])  # sort on X axis
    predict_x = [d[0] for d in data]
    predict_y = asarray([d[1] for d in data])
    plot(predict_x, predict_y, color='k', linestyle=':')