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=':')

def covariance_ci(covariance, fits, parnames, sigma=1, verbose=True):
	npars = fits.size
	out = np.zeros( (npars, 3) )
	out[:,0] = fits
	out[:,1] = fits - sigma*np.sqrt( np.diagonal(covariance) )
	out[:,2] = fits + sigma*np.sqrt( np.diagonal(covariance) )
	if sigma==1:
		percent = '68%'
	if sigma==2:
		percent = '95%'
	if verbose:
		print 'covariance '+percent+' ci'
		hfmt = '%-12s %12s %12s %12s \n'
		s = hfmt % ('Param', 'Best-Fit', 'Lower Bound', 'Upper Bound')
		s += hfmt % ('-'*5, '-'*8, '-'*11, '-'*11)
		for name, val, lower, upper in zip(parnames, out[:,0],
			out[:,1], out[:,2]):
			s += '%-12s %12g ' % (name, val)
			s += '%12g ' % lower
			s += '%12g \n' % upper
		print s
	result={}
	for i in range(npars):
		result[parnames[i]]=out[i,[1,2]]
	return result

    def calculate_energy(self):
        """
        The full energy function involves calculation of all pairwise
        energy interactions. We sum across rows down columns separately
        """
        energy = 0.
        for row in self.system:
            energy += -1 * np.sum(
                [self.coupling_constant * i * j for i,j in zip(row[:-1], row[1:])]
                )
            energy += -1 * np.sum(self.external_field * row)
            if self.next_nearest:
                energy += -1 * np.sum(			
                     [self.coupling_constant * i * j for i,j in zip(row[:-2], row[2:])]
                     )
        for col in np.rollaxis(self.system, -1):
            energy += -1 * np.sum(
                [self.coupling_constant * i * j for i,j in zip(col[:-1], col[1:])]
                )
            if self.next_nearest:
                 energy += -1 * np.sum(			
                     [self.coupling_constant * i * j for i,j in zip(row[:-2], row[2:])]
                     )
        if self.next_nearest:        #get diagnal elements 
             diags=[]
             for i in range(-len(self.system),len(self.system)):
                diags.append(np.diagonal(self.system,i,0,1))
                diags.append(np.diagonal(self.system,i,1,0))
             for diag in diags:  
                 energy += -1 * np.sum(
                    [self.coupling_constant * i * j for i,j in zip(diag[:-1], diag[1:])]
                    )

        return energy

 def __init__( 
         self, 
         data3d=None, # 3d-array,
         cell=None, # 3*3 array, vectors of cell 
         grid_measure=None, # 3*3 array, vectors of grid increament
         delta_grid=None, # 1*3 array, (dx, dy, dz)
         cell_origin=None, # 1*3 array, (x0, y0, z0)
         inptype=None, # raw input, or file type 'cube'...
         inpfname=None, # 
         ):
     '''
     Dealing with types of input:
     '''
     if inptype=='cube' and inpfname:
         self.cell_origin, self.grid_measure, self.data3d = read_cube(inpfname)
     else:
         sys.exit('Crewp:ScalarField: Input for ScalarField, not implemented.')
     '''
     Case: get ``delta_grid`` from ``grid_measure`` matrix.
     '''
     if (not hasattr(self, 'delta_grid')) and hasattr(self, 'grid_measure'):
         if np.count_nonzero(self.grid_measure - np.diag(np.diagonal(self.grid_measure))) != 0 : # test if grid_measure diagonal.
             sys.exit('Crewp:ScalarField: Non-Cartesian grid of ScalarField, not implemented.')
         else:
             self.delta_grid = np.diagonal(self.grid_measure)
     self.gridshape = self.data3d.shape

def distance_diagonal_law(matrix, positions=None):
    n = min(matrix.shape)
    if positions is None:
        return np.array([np.average(np.diagonal(matrix, j)) for j in range(n)])
    else:
        contigs = positions_to_contigs(positions)

    def is_intra(i, j):
        return contigs[i] == contigs[j]

    max_intra_distance = max((len(contigs == u) for u in set(contigs)))

    intra_contacts = []
    inter_contacts = [np.average(np.diagonal(matrix, j))
                      for j in range(max_intra_distance, n)]
    for j in range(max_intra_distance):
        D = np.diagonal(matrix, j)
        for i in range(len(D)):
            diagonal_intra = []
            if is_intra(i, j):
                diagonal_intra.append(D[i])
#            else:
#                diagonal_inter.append(D[i])
#        inter_contacts.append(np.average(np.array(diagonal_inter)))
        intra_contacts.append(np.average(np.array(diagonal_intra)))

    intra_contacts.extend(inter_contacts)

    return [positions, np.array(intra_contacts)]

 def estimate(self, X, Y, shared):
     m = len(Y)
     d = len(X.T)
     # Split spam and no spam
     nspam = np.array([x[(Y==0.0).T[0]] for x in X.T]).T
     spam = np.array([x[(Y==1.0).T[0]] for x in X.T]).T
     p = [float(len(nspam))/float(len(X)), float(len(spam))/float(len(X))]
     # Estimate u0, u1
     u0, u1 = nspam.mean(0), spam.mean(0)
     # Initialize covariance
     u = (u0, u1)
     if shared:
         cov = np.zeros((d,d))
     else:
         cov = [np.zeros((d,d)), np.zeros((d,d))]
     # Estimate covariance
     for i in range(m):
         d = (X[i] - u[int(Y[i])])
         if shared:
             cov += np.dot(d[:,None],d[None,:])
         else:
             cov[int(Y[i])] += np.dot(d[:,None],d[None,:])
     # Normalize covariance
     if shared:
         cov /= float(m)
         #print np.diagonal(cov)
         np.fill_diagonal(cov, np.diagonal(cov)+1e-5)
         #print np.diagonal(cov)
     else:
         cov = [c/float(m) for c in cov]
         [np.fill_diagonal(c, np.diagonal(c)+1e-5) for c in cov]
     # Print stats and return
     return u,cov,p