def __init__(self, mean, sigma):
        """
        @param mean: mean of Gaussian density 
        @type mean: float or array of type 'd'

        @param sigma: sigma of Gaussian density
        @type sigma: float or array of type 'd'
        """
        self.sigma=sigma
        self.mean=mean
        if type(mean)==FloatType and type(sigma)==FloatType:
            # 1D Gaussian
            self.c=-log(sigma*sqrt(2*pi))
            self.c2=2*sigma*sigma
            self.one_dim=1
        else:
            # Multidimensional Gaussian
            assert(det(sigma)>0)
            self.dim=mean.shape[0]
            assert(sigma.shape==(self.dim,self.dim))
            sq_det=sqrt(det(self.sigma))
            self.b=1.0/(power(2*pi, self.dim/2.0)*sq_det)
            self.logb=log(self.b)
            self.inv_sigma=inv(self.sigma)
            self.one_dim=0

def naiveBayes():
    global tSRiver, tSLand, im, outimg
    covRiver, covLand = cov(tSRiver, rowvar=0), cov(tSLand, rowvar=0)
    icr, icl= inv(covRiver), inv(covLand)
    meanRiver, meanLand = mean(tSRiver, axis=0), mean(tSLand, axis=0)

    for i in range(512):
        for j in range(512):
            devRiver, devLand = subtract(im[i,j], meanRiver), subtract(im[i,j], meanLand)
            rivClass = dot(dot(devRiver, icr), devRiver.T)
            nrivClass = dot(dot(devLand, icl), devLand.T)

            try:
                p1 = Decimal(-0.5)/Decimal(m.sqrt(det(covRiver))) * Decimal(m.exp(rivClass))
            except OverflowError:
                outimg[i,j] = 255
                continue# as e: print(e, 'Variable in exp: ', rivClass)
            try:
                p2 = Decimal(-0.5)/Decimal(m.sqrt(det(covLand))) * Decimal(m.exp(nrivClass))
            except OverflowError: continue# as e:
#                print(e, 'Variable in exp: ', rivClass, ' Setting pixel to 0')
            class1 = Decimal(P1)*p1
            class2 = Decimal(P2)*p2

            if class1 >= class2: outimg[i,j] = 255

def genGraph(S_actual, S_est, S_previous, empCov_set, nodeID, e1, e2, e3, e4, display = False):
    D = np.where(S_est != 0)[0].shape[0]#len(numpy.where(S_est == 0)[0])
    T = np.where(S_actual != 0)[0].shape[0]
#            print np.where(S_actual != 0)[0]
    TandD = float(np.where(np.logical_and(S_actual,S_est) == True)[0].shape[0])
    P = TandD/D
    R = TandD/T
    offDiagDiff = S_actual - S_est
    offDiagDiff = offDiagDiff - np.diag(np.diag(offDiagDiff))
    S_diff = (S_est - S_previous)  
    S_diff = S_diff - np.diag(np.diag(S_diff))
    ind = (S_diff < 1e-2) & (S_diff > - 1e-2)
    S_diff[ind] = 0    
    K = np.count_nonzero(S_diff)
    e1.append(-np.log(alg.det(S_est)) + np.trace(np.dot(S_est, empCov_set[nodeID])) + K)
#    e1.append( alg.norm(offDiagDiff, 'fro'))
    e2.append(2* P*R/(P+R))
    
    
    K = float(np.where(np.logical_and((S_est>0) != (S_previous>0), S_est>0) == True)[0].shape[0])
    e3.append(-np.log(alg.det(S_est)) + np.trace(np.dot(S_est, empCov_set[nodeID])) + K)
    e4.append(alg.norm(S_est -  S_previous, 'fro'))
    

    if display == True:
        if (nodeID >timeShift -10) and (nodeID < timeShift + 10):
            print 'nodeID = ', nodeID
            print 'S_true = ', S_actual,'\nS_est', S_est
#            print 'S_error = ',S_actual - S_est, '\n its Fro error = ', alg.norm(S_actual - S_est, 'fro')
            print 'D = ',D,'T = ', T,'TandD = ', TandD,'K = ', K,'P = ', P,'R = ', R,'Score = ', 2* P*R/(P+R)
            
    return e1, e2, e3, e4

def kabsch(mol0, mol1):
    # translate to align centroids with origin
    mol0, d0 = setPdistance(mol0, mol0.centersym(), [0, 0, 0], 0)
    mol1, d1 = setPdistance(mol1, mol1.centersym(), [0, 0, 0], 0)
    # get coordinates and matrices P,Q
    P, Q = [], []
    for atom0, atom1 in zip(mol0.getAtoms(), mol1.getAtoms()):
        P.append(atom0.coords())
        Q.append(atom1.coords())
    # Computation of the covariance matrix
    C = dot(transpose(P), Q)
    # Computation of the optimal rotation matrix
    # This can be done using singular value decomposition (SVD)
    # Getting the sign of the det(V)*(W) to decide
    # whether we need to correct our rotation matrix to ensure a
    # right-handed coordinate system.
    # And finally calculating the optimal rotation matrix U
    # see http://en.wikipedia.org/wiki/Kabsch_algorithm
    V, S, W = svd(C)
    d = (det(V) * det(W)) < 0.0
    # Create Rotation matrix U
    if d:
        S[-1] = -S[-1]
        V[:, -1] = -V[:, -1]
    U = dot(V, W)
    # Rotate P
    P = dot(P, U)
    # write back coordinates
    for i, atom in enumerate(mol0.getAtoms()):
        atom.setcoords(P[i])
    return mol0, U.tolist(), d0, d1

    def to_fsl(self, infile):
        """
        Converts a pycortex transform to an FSL transform.
        The resulting FSL transform goes FROM the space of the "infile" input
        TO the space of the reference nifti stored in the pycortex transform.

        This should ONLY be used for "coord" transforms! Will fail hard for 
        "magnet" transforms!
        """
        import nibabel
        import numpy.linalg as npl

        try:
            inIm = nibabel.load(infile)
        except AttributeError:
            inIm = infile
        in_hdr = inIm.get_header()
        ref_hdr = self.reference.get_header()
        # get_zooms gets the positive voxel sizes as returned in the header
        inspace = np.diag(in_hdr.get_zooms()[:3] + (1,))
        refspace = np.diag(ref_hdr.get_zooms()[:3] + (1,))
        # Since FSL does not use the full transform info in the nifti header, 
        # determine whether the transform indicates that the X axis should be 
        # flipped; if so, flip the X axis (for both infile and reffile)
        if npl.det(in_hdr.get_best_affine())>=0:
            print("Determinant is > 0: FLIPPING!")
            inspace = np.dot(inspace, _x_flipper(in_hdr.get_data_shape()[0]))
        if npl.det(ref_hdr.get_best_affine())>=0:
            print("Determinant is > 0: FLIPPING!")
            refspace = np.dot(refspace, _x_flipper(ref_hdr.get_data_shape()[0]))

        inAffine = inIm.get_affine()
        inv = np.linalg.inv
        fslx = np.dot(refspace,np.dot(self.xfm,np.dot(inAffine,inv(inspace))))
        return fslx

 def eval_g(self, x):
     Q = self.eval_Q(x)
     R = self.eval_R(x)
     g_Q = la.det(Q)
     g_R = la.det(R)
     g_dual = self.eval_dual_obj(x)
     return np.array([g_Q, g_R, g_dual])

    def from_fsl(cls, xfm, infile, reffile):
        """
        Converts fsl transform xfm (estimated FROM infile TO reffile) 
        to a pycortex COORD transform. 
        """
        ## Adapted from dipy.external.fsl.flirt2aff#############################
        import nibabel
        import numpy.linalg as npl
        
        try:
            inIm = nibabel.load(infile)
        except AttributeError:
            inIm = infile
        
        refIm = nibabel.load(reffile)
        in_hdr = inIm.get_header()
        ref_hdr = refIm.get_header()
        # get_zooms gets the positive voxel sizes as returned in the header
        inspace = np.diag(in_hdr.get_zooms()[:3] + (1,))
        refspace = np.diag(ref_hdr.get_zooms()[:3] + (1,))
        # Since FSL does not use the full transform info in the nifti header, 
        # determine whether the transform indicates that the X axis should be 
        # flipped; if so, flip the X axis (for both infile and reffile)
        if npl.det(in_hdr.get_best_affine())>=0:
            inspace = np.dot(inspace, _x_flipper(in_hdr.get_data_shape()[0]))
        if npl.det(ref_hdr.get_best_affine())>=0:
            refspace = np.dot(refspace, _x_flipper(ref_hdr.get_data_shape()[0]))

        inAffine = inIm.get_affine()
        inv = np.linalg.inv
        coord = np.dot(inv(refspace),np.dot(xfm,np.dot(inspace,inv(inAffine))))
        return cls(coord, refIm)

 def marginal(self, X, a, B, c, m):
     d = float(X.shape[0])
     px1 = (c / (math.pi * (1 + c)))**(d/2)
     px2 = det(B + c/(1+c) * (X - m) * (X - m).T)**(-(a+1)/2) / det(B)**(-a/2)
     px3 = math.exp(np.sum( math.log(gamma( (a+1)/2 + (1 - np.arange(1.0, d))/2 )) - math.log(gamma( a/2 + (1 - np.arange(1.0, d))/2 )) ))
     px = px1 * px2 * px3
     return px

def diwishart (X,v,S,return_log=False):
    is_square(X)
    k = X.ndim
    if not return_log: 
        return np.exp(-0.5*(S*la.inv(X)).trace()) * pow(la.det(X),(-(v+k+1)/2.0)) * pow(la.det(S),(v/2.0)) / ( pow(2,(v*k/2.0)) * pow(pi,(k*(k-1)/4.0))* gamma((v+1-np.array(range(1,k)))/2.0).prod())
    else:
        return -0.5*(S*la.inv(X)).trace() + np.log(la.det(X))*(-(v+k+1)/2.) +np.log(la.det(S))*(v/2.0) - np.log(2)*(v*k/2.0) -np.log(pi)*(k*(k-1)/4.0) - lgamma((v+1-np.array(range(1,k)))/2.0).sum()

def findmin(M,B,run): #normal routine for varying a single element at a time. 
    '''Finds minimum cost for the lattice by varying the integers of m, in the element that gives steepest descent
    The 'run' indicates the cost function to use'''
    maxsteps = 10000
    istep = 1
    M = minkM(M,B)#;print'Mink reduced M'; print M
    
    while istep<maxsteps:

#        sys.exit('stop')
        bestindex = changewhich(M,B,run)
        if bestindex[2]==0:#found minimum
#            newcost = cost(M,B,run)
            break
        else:
            M[bestindex[0],bestindex[1]] += bestindex[2]
#            newcost = cost(M,B,run)
            if run == 'minsvsym' and B.Nmesh/float(det(M))>1.2: M = M*2 # open up search space when when det(M) gets too low
#            print; print M;print newcost

        istep += 1
    if istep < maxsteps:
#        print 'Found minimum after %i steps' % istep
#        print 'Best M'; print M
        K = lattice();K.vecs = trimSmall(dot(B.vecs,inv(M)));K.det = abs(det(K.vecs)); K.Nmesh = B.det/K.det             
    else:
        print 'Ended without minimum after maximum %i steps' % istep
        sys.exit('Stop')
    return [trimSmall(M),K]

def nk_bhatta(X1, X2, eta):
    # Returns the non-kernelized Bhattacharrya
    #I.e. fits normal distributions in input space and calculates Bhattacharrya overlap between them
    (n1, d1) = X1.shape
    (n2, d ) = X2.shape
    assert d1 == d
    mu1 = mat.sum(X1,0) / n1
    mu2 = mat.sum(X2,0) / n2
    X1c = X1 - mat.tile(mu1, (n1,1))
    X2c = X2 - mat.tile(mu2, (n2,1))
    Eta = mat.eye(d) * eta
    S1 = X1c.T * X1c / n1 + Eta
    S2 = X2c.T * X2c / n2 + Eta

    mu3 = .5 * (S1.I * mu1.T + S2.I * mu2.T).T
    S3  = 2  * (S1.I + S2.I).I

    d1 = la.det(S1) ** -.25
    d2 = la.det(S2) ** -.25
    d3 = la.det(S3) ** .5
    dterm = d1 * d2 * d3

    e1 = -.25 * mu1 * S1.I * mu1.T
    e2 = -.25 * mu2 * S2.I * mu2.T
    e3 = .5   * mu3 * S3   * mu3.T

    eterm = math.exp(e1 + e2 + e3)

    return float(dterm * eterm)

def plot_Peierls_Hubbard_green(u=0.0):
	''' plot the Det(G) in the plane dt-w'''
	L = 6
	num = 51
	d_num = 20
	ws = np.linspace(-3,3,num)
	Gd = np.zeros((num, d_num), dtype=np.complex)
	for j, d in enumerate(np.linspace(-1,1,d_num)):
		H = Peierls_Hubbard(L, t=1., dt=d, U=u, pbc=False)
		evals, evecs = func.solve_hamiltonian(H)

		Gwk = np.zeros(num, dtype=np.complex)
		Gij_up = func.green(ws, op.fermion_up(), L, evals, evecs)
		Gij_down = func.green(ws, op.fermion_down(), L, evals, evecs)
		for i in range(num):
			Gwk[i] =  det(Gij_up[i])*det(Gij_down[i])
		Gd[:,j] = Gwk

		# plt.figure("Green u=%.2f"%d)
		# plt.plot(ws, np.real(Gwk),'o-')
		# plt.ylim([-1,1])
		# print np.min( np.imag(Gwk) ), np.max( np.imag(Gwk) )

	plt.figure("poles")
	fig = plt.imshow(np.real(Gd), interpolation='nearest', aspect=0.4)
	fig.set_clim(-1,1)
	plt.colorbar(ticks=[-1,0,1],label="Det(G)")
	plt.xticks([0.5,d_num/2.,d_num-0.5], [-1,0,1], fontsize=20)
	plt.yticks([0.5,num/2.,num-0.5], [3,0,-3], fontsize=20)
	plt.xlabel('d', fontsize=20)
	plt.ylabel('w', fontsize=20)

 def update(self, x, y):
     tmp_X = self.X
     tmp_y = self.y
     tmp_reg_K = self.reg_K
     
     self.X = vstack((self.X, x))
     self.y = hstack((self.y, y))
     
     n = self.X.shape[0]
     if (n <= th):
         if (n != 1):
             row = zeros((1, n-1))
             col = zeros((n, 1))
             self.reg_K = vstack((self.reg_K, row))
             self.reg_K = hstack((self.reg_K, col))
         self.add_b(x, n)
         if det(self.reg_K) == 0:
             self.X = tmp_X
             self.y = tmp_y
             self.reg_K = tmp_reg_K
             return
         self.reg_K_inv = self.reg_K.getI()
     else:
         self.X = delete(self.X, 0, 0)
         self.y = delete(self.y, 0, 0)
         self.updateK(x)
         if det(self.reg_K) is 0:
             self.X = tmp_X
             self.y = tmp_y
             self.reg_K = temp_reg_K
             return
         self.updateKInv()

def circleFromPoints(p1, p2, p3):
    """
    This function takes three points in the image and returns the center and radius of the circle passing through these 
    three points
    """
    x1, y1 = p1
    x2, y2 = p2
    x3, y3 = p3
    centerx = det(np.array(([x1**2+y1**2,y1,1],[x2**2+y2**2,y2,1],[x3**2+y3**2,y3,1]))) / float(2*det(np.array(([x1,y1,1],[x2,y2,1],[x3,y3,1]))))
    centery = det(np.array(([x1,x1**2+y1**2,1],[x2,x2**2+y2**2,1],[x3,x3**2+y3**2,1]))) / float(2*det(np.array(([x1,y1,1],[x2,y2,1],[x3,y3,1]))))
    radius = sqrt( pow(x2 - centerx, 2) + pow(y2-centery, 2))
    # Alternative:
    # offset = pow(x2, 2) +pow(y2, 2)
    # bc = ( pow(x1, 2) + pow(y1, 2) - offset )/2.0
    # cd = (offset - pow(x3, 2) - pow(y3, 2))/2.0
    # det = (x1 - x2) * (y2 - y3) - (x2 - x3)* (y1 - y2) 
    # TOL = 0.0000001
    # if (abs(det) < TOL):
    #     print "POINTS TOO CLOSE"
    #     return ((0, 0), 0) 
    # idet = 1/det
    # centerx =  (bc * (y2 - y3) - cd * (y1 - y2)) * idet
    # centery =  (cd * (x1 - x2) - bc * (x2 - x3)) * idet
    # print centerx, centery, radius
    return ((centerx, centery), radius)

    def InitializeElectrons(self):
        if self.N_up == 0 and self.N_down == 0:
            print 'Error: no electrons to initialize'
            return []
        else:
            # generate array of electron positions, normally distributed from the origin with Bohr radius
            n = self.N_up+self.N_down
            self.e_positions = np.random.randn(n,3) * GSF.a_B # generate array of electron positions
            print 'init e_pos',self.e_positions
            # Store displacements and distances
            self.e_disp = np.zeros((n,n,3)) # store the displacements in a 3D matrix to make indexing easier
            self.e_dist = np.zeros((n,n)) # the electron matrices should only be upper diagonal
            self.atom_disp = np.zeros((n,len(self.atom_list),3))
            self.atom_dist = np.zeros((n,len(self.atom_list)))
            index = 0
            for i in range(n):
                self.e_disp[i,i+1:] = self.e_positions[i] - self.e_positions[i+1:]
                self.e_dist[i,i+1:] = np.sqrt(np.sum(self.e_disp[i,i+1:]**2,1))
                self.atom_disp[i] = self.e_positions[i] - self.ion_positions
                self.atom_dist[i,:] = np.sqrt(np.sum(self.atom_disp[i,:]**2,1))
 	      #Once the e_position is initialize, the slater matrix and its deteriminant and inverse are all initialized. 
        self.slater_matrix_up = SlaterMatrix(self.e_positions[0:self.N_up],self.psi_up)
        self.slater_matrix_down = SlaterMatrix(self.e_positions[self.N_up:],self.psi_down)
        print 'slater_matrix',self.slater_matrix_up
        if self.N_up>0: 
            self.inverse_SD_up = LA.inv(self.slater_matrix_up)
            self.slater_det_up = LA.det(self.slater_matrix_up)			
        if self.N_down>0: 
            self.inverse_SD_down = LA.inv(self.slater_matrix_down) 
            self.slater_det_down = LA.det(self.slater_matrix_down)
        self.J = self.Jastrow()
        print 'slater_inv',self.inverse_SD_up
        return self.e_positions

def marginal_log_likelihood(y, A,Q,C,R, _x,_P):
    """
    compute the marginal log likelihood
    given the learned model parameters
    and the learned expected sufficient statistics

    C must be invertible (:= square and non-singular)
    eg aint no inverse to a fucking projection

    p(y|params) = ...
    """

    if True: return 0

    T,p = y.shape
    _,d = _x.shape
    Cy = mul(inv(C), tp(y))
    CRC = inv(mul(tp(C),inv(R),C))

    Z_swap = (1/2) * (log(det(2*pi*CRC)) - log(det(2*pi*R)))

    covar = CRC+identity(d)
    Z_merge = log(2*pi*det(covar))**(-d/2) - (1/2)*mul(tp(Cy), covar, Cy)

    for t in r(1,T-1):
        mean  = _x[t+1]
        covar = CRC + _P[t+1]
        Z_merge += log(2*pi*det(covar))**(-d/2) - (1/2)*mul(tp(Cy-mean), covar, (Cy-mean))

    return sum(T*Z_swap + Z_merge)

 def _calculate_log_likelihood(self):
     #if self.m == None:
     #    Give error message
     R = zeros((self.n, self.n))
     X,Y = array(self.X), array(self.Y)
     thetas = 10.**self.thetas
     for i in range(self.n):
         for j in arange(i+1,self.n):
             R[i,j] = (1-self.nugget)*e**(-sum(thetas*(X[i]-X[j])**2.)) #weighted distance formula
     R = R + R.T + eye(self.n)
     self.R = R
     one = ones(self.n)
     try:
         self.R_fact = cho_factor(R)
         rhs = vstack([Y, one]).T
         R_fact = (self.R_fact[0].T,not self.R_fact[1])
         cho = cho_solve(R_fact, rhs).T
         
         self.mu = dot(one,cho[0])/dot(one,cho[1])
         self.sig2 = dot(Y-dot(one,self.mu),cho_solve(self.R_fact,(Y-dot(one,self.mu))))/self.n
         #self.log_likelihood = -self.n/2.*log(self.sig2)-1./2.*log(abs(det(self.R)+1.e-16))-sum(thetas)
         self.log_likelihood = -self.n/2.*log(self.sig2)-1./2.*log(abs(det(self.R)+1.e-16))
     except (linalg.LinAlgError,ValueError):
         #------LSTSQ---------
         self.R_fact = None #reset this to none, so we know not to use cholesky
         #self.R = self.R+diag([10e-6]*self.n) #improve conditioning[Booker et al., 1999]
         rhs = vstack([Y, one]).T
         lsq = lstsq(self.R.T,rhs)[0].T
         self.mu = dot(one,lsq[0])/dot(one,lsq[1])
         self.sig2 = dot(Y-dot(one,self.mu),lstsq(self.R,Y-dot(one,self.mu))[0])/self.n
         self.log_likelihood = -self.n/2.*log(self.sig2)-1./2.*log(abs(det(self.R)+1.e-16))

def compDklgaussian(mu1, C1, mu2, C2):
    '''
    mu1 and C1 specify a multidimensional gaussian. mu2 and C2 specify another
    one (of the same dimension). Return is a float giving the KL divergence
    between the two Gaussians.

    Inputs have to be matrix instances, and the mu inputs should be in row
    vector shape (mu.shape[0] = 1, mu.shape[1] > 1)

    '''
    n = mu1.size
    b = mu2 - mu1
    C2inv = la.inv(C2)
    C1sqrt = np.mat(sqrtPSDm(C1))
    Term1 = C1sqrt * C2inv * C1sqrt
    Term2 = b.transpose() * C2inv * b
    det1 = la.det(C1)
    det2 = la.det(C2)
    tol = 1e8
    if (la.cond(C1) > tol) | (la.cond(C2) > tol):
        print('Determinants not non-zero. Ignoring determinants.')
        Term3 = 0
    else:
        Term3 = .5 * np.log(det2 / det1)
    d = .5 * np.trace(Term1) + .5 * Term2 - .5 * n + Term3
    return d[0, 0]

def read_data_inhomogeneous(iom, blockid=0):
    r"""
    :param iom: An :py:class:`IOManager` instance providing the simulation data.
    :param blockid: The data block from which the values are read.
    """
    parameters = iom.load_parameters()
    timegrid = iom.load_inhomogwavepacket_timegrid(blockid=blockid)
    dt = parameters["dt"] if "dt" in parameters else 1.0
    # Filter
    time = timegrid * dt
    time = where(timegrid < 0, nan, time)

    Pis = iom.load_inhomogwavepacket_parameters(blockid=blockid)

    # List with Pi time evolutions
    Phist = []
    Qhist = []
    Shist = []
    phist = []
    qhist = []

    for q, p, Q, P, S in Pis:
        qhist.append(array([norm(q[i, :]) for i in range(q.shape[0])]).reshape(-1))
        phist.append(array([norm(p[i, :]) for i in range(p.shape[0])]).reshape(-1))
        Qhist.append(array([det(Q[i, :, :]) for i in range(Q.shape[0])]).reshape(-1))
        Phist.append(array([det(P[i, :, :]) for i in range(P.shape[0])]).reshape(-1))
        Shist.append(S.reshape(-1))

    return (time, qhist, phist, Qhist, Phist, Shist)

def ComputeIndexLda(proj, data, dataLabels): 
# takken from page 31 of book on GGOBI
# diane cook "integrative and dynamic graphics for data analysis
    ndim = proj.shape[0]

    A = matrix(proj)
    B = matrix(zeros((ndim, ndim)))
    W = matrix(zeros((ndim, ndim)))
    grand_mean = mean(data , axis = 0).astype(matrix)
    for c in unique(dataLabels):
        ind = dataLabels ==c
        group_mean = mean(data[ind, :], axis = 0).astype(matrix)
        d = group_mean - grand_mean
        B += ind[ind].size*dot(d[:,newaxis], d[newaxis, :])
        
        # find better:
        d = data[ind, :] - group_mean
        for j in range(data[ind, :].shape[0]):
            W+=dot(d[j,:, newaxis], d[j,newaxis,:])
        
    #~ print 'A', A, A.size
    #~ print 'B', B, B.size
    #~ print 'W', W, W.size

    lda = 1 -  linalg.det(A.T*W*A) / linalg.det(A.T*(W+B)*A)
    return lda