def berrycurvature(H,kx,ky,mu,d=10**-6):
    '''
    Calculate the Berry curvature of the occupied bands for a Hamiltonian with the given chemical potential using the Kubo formula.

    Parameters
    ----------
    H : callable
        Input function which returns the Hamiltonian as a 2D array.
    kx,ky : float
        The two parameters which specify the 2D point at which the Berry curvature is to be calculated.
        They are also the input parameters to be conveyed to the function H.
    mu : float
        The chemical potential.
    d : float, optional
        The spacing to be used to calculate the derivatives.

    Returns
    -------
    float
        The calculated Berry curvature for function H at point kx,ky with chemical potential mu.
    '''
    result=0
    Vx=(H(kx+d,ky)-H(kx-d,ky))/(2*d)
    Vy=(H(kx,ky+d)-H(kx,ky-d))/(2*d)
    Es,Evs=eigh(H(kx,ky))
    for n in range(Es.shape[0]):
        for m in range(Es.shape[0]):
            if Es[n]<=mu and Es[m]>mu:
                result-=2*(np.vdot(np.dot(Vx,Evs[:,n]),Evs[:,m])*np.vdot(Evs[:,m],np.dot(Vy,Evs[:,n]))/(Es[n]-Es[m])**2).imag
    return result

    def testMCMCFitting(self):
        print "starting mcmc"
        true, data = self.generate_data()
        model = HDPMixtureModel(3,100,100,1)
        model.seed = 1
        start = time()
        r = model.fit(data, verbose=10)
#        end = time() - start
#        print r.pis.shape
#        print r.pis, r.pis[0]
#        print r[0],r[0].pis
        print r.mus
        self.assertEqual(len(r), 2, 'results are the wrong length: %d' %len(r))
        diffs = {}
        #print 'r.mus:', r.mus()
        for i in gen_mean:
            diffs[i] = np.min(np.abs(r[0].mus-gen_mean[i]),0)
            #print r[0].mus[0], np.min(np.abs(r[0].mus[0]-gen_mean[i]),0)
            #diffs[i] = np.abs(r[0].mus()[i]-gen_mean[i])
            #print i, gen_mean[i],r[0].mus()[i], diffs[i], np.vdot(diffs[i],diffs[i])
            self.assertLessEqual( np.vdot(diffs[i],diffs[i]),1, 
                                  'diff to large: %f' % np.vdot(diffs[i], diffs[i]))
        
        for i in gen_mean:
            diffs[i] = np.min(np.abs(r[1].mus-gen_mean[i]),0)
            #diffs[i] = np.abs(r[1].mus()[i]-gen_mean[i])
            #print i, gen_mean[i],r[1].mus()[i], diffs[i], np.vdot(diffs[i],diffs[i])
            self.assertLessEqual( np.vdot(diffs[i],diffs[i]),1, 
                                  'diff to large: %f' % np.vdot(diffs[i], diffs[i]))

    def test_integral_moment_first(self):
        lmax = 5 #XXX use meshgrid above l=5

        # Test first-order moments of the spherical harmonics
        for l1,m1 in lmiter(lmax, comm=world):
            s1_L = Y(l1, m1, theta_L, phi_L)
            for l2,m2 in lmiter(lmax):
                s2_L = Y(l2, m2, theta_L, phi_L)

                # Note that weights times surface area make up for sin(theta)
                v_ex = 4 * np.pi * np.vdot(s1_L, np.cos(phi_L) \
                    * np.sin(theta_L) * s2_L * weight_L)
                v_ey = 4 * np.pi * np.vdot(s1_L, np.sin(phi_L) \
                    * np.sin(theta_L) * s2_L * weight_L)
                v_ez = 4 * np.pi * np.vdot(s1_L, np.cos(theta_L) \
                    * s2_L * weight_L)

                v0_ex = intYY_ex(l1, m1, l2, m2)
                v0_ey = intYY_ey(l1, m1, l2, m2)
                v0_ez = intYY_ez(l1, m1, l2, m2)

                self.assertAlmostEqual(v_ex, v0_ex, 12, '%s != %s (l1=%2d, ' \
                    'm1=%2d, l2=%2d, m2=%2d)' % (v_ex,v0_ex,l1,m1,l2,m2))
                self.assertAlmostEqual(v_ey, v0_ey, 12, '%s != %s (l1=%2d, ' \
                    'm1=%2d, l2=%2d, m2=%2d)' % (v_ey,v0_ey,l1,m1,l2,m2))
                self.assertAlmostEqual(v_ez, v0_ez, 12, '%s != %s (l1=%2d, ' \
                    'm1=%2d, l2=%2d, m2=%2d)' % (v_ez,v0_ez,l1,m1,l2,m2))

def neb_modify_force(band):
	peak_image = numpy.argmax(band.Energy)
	for i in xrange(1, band.nimages-1):
		# "improved" tangent calculation
		Forwards = putil.separation_vector_raw(band.Pos[(i+1)*3*band.NAtoms:(i+2)*3*band.NAtoms], band.Pos[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], band.PBC, band.Dim)
		Backwards = putil.separation_vector_raw(band.Pos[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], band.Pos[(i-1)*3*band.NAtoms:(i)*3*band.NAtoms], band.PBC, band.Dim)
		mForwards = putil.magnitude(Forwards)
		mBackwards = putil.magnitude(Backwards)
	
		if(band.Energy[i]>band.Energy[i-1] and band.Energy[i+1] > band.Energy[i]):
			norm = putil.normalise(Forwards)
		elif(band.Energy[i] < band.Energy[i-1] and band.Energy[i+1] < band.Energy[i]):
			norm = putil.normalise(Backwards)
		else:
			norm = putil.normalise(Backwards * math.fabs(band.Energy[i] - band.Energy[i-1]) + Forwards * math.fabs(band.Energy[i+1] - band.Energy[i]))

			
		# if climbing then move uphill independent of the rest of the band
		if(i==peak_image):
			band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms] -= 2 * norm * numpy.vdot(band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], norm)
			
		else:
			# remove parallel component of real force
			band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms] -= norm * numpy.vdot(band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], norm)
			# add in force due to the springs
			band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms] += norm*(mForwards - mBackwards)*band.springk
			
		#band.mForce = numpy.absolute(band.Force).max()
		band.mForce = putil.magnitude(band.Force)
		band.peak_image = peak_image
		

		
	return band	

def adjointness_error(op, its=100):
    """Check adjointness of op.A and op.As for 'its' instances of random data.

    For random unit-normed x and y, this finds the error in the adjoint
    identity <Ax, y> == <x, A*y>:
        err = abs( vdot(A(x), y) - vdot(x, Astar(y)) ).

    The type and shape of the input to A are specified by inshape and indtype.

    Returns a vector of the error magnitudes.

    """
    inshape = op.inshape
    indtype = op.indtype
    outshape = op.outshape
    outdtype = op.outdtype

    x = get_random_normal(inshape, indtype)

    errs = np.zeros(its, dtype=indtype)
    for k in xrange(its):
        x = get_random_normal(inshape, indtype)
        x = x/np.linalg.norm(x)
        y = get_random_normal(outshape, outdtype)
        y = y/np.linalg.norm(y)
        ip_A = np.vdot(op.A(x), y)
        ip_Astar = np.vdot(x, op.As(y))
        errs[k] = np.abs(ip_A - ip_Astar)

    return errs

    def braket(self, other, space="g"):
        """
        Returns the scalar product <u1|u2> of the periodic part of two wavefunctions 
        computed in G-space or r-space, depending on the value of space.

        Args:
            other: 
                Other wave (right-hand side)
            space: 
                Integration space. Possible values ["g", "gsphere", "r"]
                if "g" or "r" the scalar product is computed in G- or R-space on the FFT box.
                if space="gsphere" the integration is done on the G-sphere. Note that
                this option assumes that self and other have the same list of G-vectors. 
        """
        space = space.lower()

        if space == "g":  
            ug1_mesh = self.gsphere.tofftmesh(self.mesh, self.ug)
            ug2_mesh = other.gsphere.tofftmesh(self.mesh, other.ug)
            return np.vdot(ug1_mesh, ug2_mesh)

        elif space == "gsphere":
            return np.vdot(self.ug, other.ug)

        elif space == "r":
            return np.vdot(self.ur, other.ur) * self.mesh.dv

        else:
            raise ValueError("Wrong space: %s" % space)

    def  __compute_velocities_at_collision(obj1, obj2):

        v_obj1 = np.array([obj1.velocity.x, obj1.velocity.y], dtype=np.float64)
        v_obj2 = np.array([obj2.velocity.x, obj2.velocity.y], dtype=np.float64)

        v_normal = np.array([v_obj1[0] - v_obj2[0], v_obj1[1] - v_obj2[1]])
        #print v_obj1, v_obj2
        #print v_normal

        if not np.linalg.norm(v_normal):
            #print "No colision..."
            return obj1.velocity, obj2.velocity

        v_unit_normal = np.divide(v_normal, np.linalg.norm(v_normal))
        v_unit_tangent = np.array([-1*v_unit_normal[1], v_unit_normal[0]])

        v_obj1_normal = np.vdot(v_unit_normal, v_obj1)
        v_obj1_tangent = np.vdot(v_unit_tangent, v_obj1)
        v_obj2_normal = np.vdot(v_unit_normal, v_obj2)
        v_obj2_tangent = np.vdot(v_unit_tangent, v_obj2)

        v_obj1_new_normal = (v_obj1_normal*(obj1.mass - obj2.mass) + 2*obj2.mass*v_obj2_normal)/(obj1.mass + obj2.mass)
        v_obj2_new_normal = (v_obj2_normal*(obj2.mass - obj1.mass) + 2*obj1.mass*v_obj1_normal)/(obj1.mass + obj2.mass)

        vec_obj1_normal = np.dot(v_obj1_new_normal, v_unit_normal)
        vec_obj1_tangent = np.dot(v_obj1_tangent, v_unit_tangent)
        vec_obj2_normal = np.dot(v_obj2_new_normal, v_unit_normal)
        vec_obj2_tangent = np.dot(v_obj2_tangent, v_unit_tangent)

        v_obj1_final = np.add(vec_obj1_normal, vec_obj1_tangent)
        v_obj2_final = np.add(vec_obj2_normal, vec_obj2_tangent)

        return Velocity(v_obj1_final[0], v_obj1_final[1]), Velocity(v_obj2_final[0], v_obj2_final[1])

    def step(self,r,f):
        r = np.array(r)
        f = np.reshape(f,(-1,3))
        if self.v is None:
            self.v = np.zeros((len(f), 3))
        else:
            try:
                vf = np.vdot(self.v,f)
            except ValueError:
                self.v = np.zeros((len(f), 3))
                vf = np.vdot(self.v,f)

            if vf > 0.0:
                self.v = ((1.0 - self.a) * self.v + 
                          (self.a * f / 
                           np.sqrt(np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v))))
                if self.Nsteps > self.Nmin:
                    self.dt = min(self.dt * self.finc, self.dtmax)
                    self.a *= self.fa
                self.Nsteps += 1
            else:
                self.v *= 0.0
                self.a = self.astart
                self.dt *= self.fdec
                self.Nsteps = 0

        self.v += self.dt * f
        dr = self.dt * self.v
        dr = self.determine_step(dr)
        return r + dr.reshape(r.shape)

def plot_ritz(A, n, iters):
    ''' Plot the relative error of the Ritz values of `A'.
    '''
    Amul = A.dot
    b = np.random.rand(A.shape[0])
    Q = np.empty((len(b), iters+1), dtype = np.complex128)
    H = np.zeros((iters+1, iters), dtype = np.complex128)
    Q[:, 0] = b / la.norm(b)
    eigvals = np.sort(abs(la.eig(A)[0]))[::-1]
    eigvals = eigvals[:n]
    abs_err = np.zeros((iters,n))

    for j in xrange(iters):
        Q[:, j+1] = Amul(Q[:, j])
        for i in xrange(j+1):
            H[i,j] = np.vdot(Q[:,i].conjugate(), (Q[:, j+1]))
            Q[:,j+1] = Q[:,j+1] - H[i,j] * (Q[:,i])

        H[j+1, j] = np.sqrt(np.vdot(Q[:, j+1], Q[:, j+1].conjugate()))
        Q[:,j+1] = Q[:,j+1] / H[j+1, j]

        if j < n:
            rit = np.zeros(n, dtype = np.complex128)
            rit[:j+1] = np.sort(la.eig(H[:j+1, :j+1])[0])[::-1]
            abs_err[j,:] = abs(eigvals - rit) / abs(eigvals)
        else:
            rit = np.sort(la.eig(H[:j+1,:j+1])[0])[::-1]
            rit = rit[:n]
            abs_err[j,:] = abs(eigvals - rit) / abs(eigvals)

    for i in xrange(n):
        plt.semilogy(abs_err[:,i])
    plt.show()

def bilinear_concentric_potential(r_g, dr_g, f_g, ft_g, l1, l2, alpha, rfilt=None):
    """Calculate corrections for concentric functions and potentials::

                 /     _   _a    _   _a    ~   _   _a  ~ _   _a   _
        v      = | f  (r - R ) V(r - R ) - f  (r - R ) V(r - R ) dr
         m1,m2   /  L1,L2                    L1,L2

    where f(r) and ft(r) are bilinear product of two localized functions which
    are radial splines times real spherical harmonics (l1,m1) or (l2,m2) and::

          _       1       _ -1              ~ _    erf(alpha*r)  _ -1
        V(r) = --------- |r|        ^       V(r) = ------------ |r|
               4 pi eps0                            4 pi eps0

    Note that alpha (and rfilt) should conform with the cutoff radius.
    """
    work_g = erf(alpha*r_g)

    if rfilt is None:
        M = np.vdot(f_g - ft_g * work_g, r_g * dr_g)
    else:
        M = np.vdot((f_g - ft_g * work_g)[r_g>=rfilt], \
            (r_g * dr_g)[r_g>=rfilt])

        # Replace 1/r -> (3-r^2/rfilt^2)/(2*rfilt) for r < rfilt
        M += np.vdot((f_g - ft_g * work_g)[r_g<rfilt], \
            (r_g**2/(2*rfilt) * (3-(r_g/rfilt)**2) * dr_g)[r_g<rfilt])

    v_mm = np.empty((2*l1+1, 2*l2+1), dtype=float)
    for m1 in range(2*l1+1):
        for m2 in range(2*l2+1):
            v_mm[m1,m2] = M * intYY(l1, m1-l1, l2, m2-l2)
    return v_mm

    def calc_orbital_params(self):
        #Calculates Keplerian orbital parameters based on the state vectors
        #This method should be called after each change to velocity vector

        #a = - mu / (v^2 - 2 * mu/r)
        v_2 = numpy.vdot(self.v, self.v)
        r_ = numpy.linalg.norm(self.r)
        # TODO: this won't work for parabolic trajectories (as the denominator is 0)!
        self.a = -EARTH_MU / (v_2 - 2 * EARTH_MU / r_)

        #T = 2*Pi*sqrt(a^3/ni)
        # TODO: again, orbital period is not defined for non-bound orbits
        self.T = 2.0 * numpy.pi * numpy.sqrt(self.a**3 / EARTH_MU)

        #Calculate specific relative angular momentum h = r X v
        h = numpy.cross(self.r, self.v)
        h_2 = numpy.vdot(h, h)
        h_ = numpy.sqrt(h_2)

        #Calculate eccentricity vector e = (v X h) / EARTH_MU - r/|r|
        e = numpy.cross(self.v, h) / EARTH_MU - self.r/r_
        self.e = numpy.linalg.norm(e)

        i_rad = numpy.arccos(h[2] / h_) #h[2] = hz
        self.i = numpy.degrees(i_rad)
        #However, some soruces state that if hz < 0 then inclination is 180 deg - i; should check this
        #n is the vector pointing to the ascending node
        n = numpy.array((-h[1], h[0], 0))
        n_ = numpy.linalg.norm(n)
        if i_rad == 0.0:
            o_rad = 0.0
        else:
            if n[1] >= 0.0: #ie. if h[0] >= 0
                o_rad = numpy.arccos(n[0] / n_)
            else:
                o_rad = 2 * numpy.pi - numpy.arccos(n[0] / n_)
        self.o = numpy.degrees(o_rad)

        #Calculate ni (true anomaly)
        q = numpy.vdot(self.r, self.v) #What the hell is q?
        ni_x = h_2 / (r_ * EARTH_MU) - 1.0
        ni_y = h_ * q / (r_ * EARTH_MU)
        self.ni = numpy.degrees(numpy.arctan2(ni_y, ni_x))

        if self.e == 0.0:
            #For circular orbit w is 0 by convention
            self.w = 0.0
        else:
            if n_ == 0.0:
                #For equatorial orbit
                self.w = numpy.degrees(numpy.arctan2(e[1], e[0]))
            else:
                self.w = numpy.degrees(numpy.arccos(numpy.vdot(n, e) / (n_ * self.e)))
        if e[2] < 0.0:
            self.w = 360.0 - self.w
        if self.w < 0.0:
            self.w = 360.0 + self.w

        self.rp = self.a * (1.0 - self.e) #Periapsis distance
        self.ra = self.a * (1.0 + self.e) #Apoapsis distance

def CGSolve(u0x,u0y,lamb,mu,b,epsilon,dfx,dfy) :
    # Solves JTJ[ux,uy]=b
    #lambd,mu,epsilon,dfx,dfy are needed in the computation of JTJ
    # [u0x,u0y] is the starting point of the iteration algo
    nitmax=100;
    ux=u0x;
    uy=u0y;
    # Computes JTJu
    Ax,Ay=JTJ(ux,uy,dfx,dfy,lamb,mu,epsilon);
    rx=b[0]-Ax
    ry=b[1]-Ay
    px=rx
    py=ry
    rsold=np.linalg.norm(rx)**2+np.linalg.norm(ry)**2
    for i in range(nitmax) :
        Apx,Apy=JTJ(px,py,dfx,dfy,lamb,mu,epsilon);
        alpha=rsold/(np.vdot(rx[:],Apx[:])+np.vdot(ry[:],Apy[:]))
        ux=ux+alpha*px
        uy=uy+alpha*py
        rx=rx-alpha*Apx
        ry=ry-alpha*Apy
        rsnew=np.linalg.norm(rx)**2+np.linalg.norm(ry)**2
        if np.sqrt(rsnew)<1e-10 :
            return [ux,uy]
        px=rx+rsnew/rsold*px
        py=ry+rsnew/rsold*py
        rsold=rsnew
    return [ux,uy]

def berryphase(H,path,ns):
    '''
    Calculate the Berry phase of some bands of a Hamiltonian along a certain path.

    Parameters
    ----------
    H : callable
        Input function which returns the Hamiltonian as a 2D array.
    path : iterable of dict
        The path along which to calculate the Berry phase.
    ns : iterable of int
        The sequences of bands whose Berry phases are wanted.

    Returns
    -------
    1d ndarray
        The wanted Berry phase of the selected bands.
    '''
    ns=np.array(ns)
    for i,parameters in enumerate(path):
        new=eigh(H(**parameters))[1][:,ns]
        if i==0:
            result=np.ones(len(ns),new.dtype)
            evs,old=new,new
        else:
            for j in range(len(ns)):
                result[j]*=np.vdot(old[:,j],new[:,j])
            old=new
    else:
        for j in range(len(ns)):
            result[j]*=np.vdot(old[:,j],evs[:,j])
    return np.angle(result)/np.pi

    def step(self, f):
        coords = self.coords
        if self.v is None:
            self.v = np.zeros((len(coords)))
        else:
            vf = np.vdot(f, self.v)
            if vf > 0.0:
                self.v = (1.0 - self.a) * self.v + self.a * f / np.sqrt(
                    np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v))
                if self.Nsteps > self.Nmin:
                    self.dt = min(self.dt * self.finc, self.dtmax)
                    self.a *= self.fa
                self.Nsteps += 1
            else:
                self.v[:] *= 0.0
                self.a = self.astart
                self.dt *= self.fdec
                self.Nsteps = 0

        self.v += self.dt * f
        dr = self.dt * self.v
        if False:  # how do we determine maxstep?
            normdr = np.sqrt(np.vdot(dr, dr))
        else:
            normdr = max(np.abs(dr))
        if normdr > self.maxstep:
            dr = self.maxstep * dr / normdr
        self.coords = coords + dr

    def step(self,f):
        atoms = self.atoms
        if self.v is None:
            self.v = np.zeros((len(atoms), 3))
        else:
            vf = np.vdot(f, self.v)
            if vf > 0.0:
                self.v = (1.0 - self.a) * self.v + self.a * f / np.sqrt(
                    np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v))
                if self.Nsteps > self.Nmin:
                    self.dt = min(self.dt * self.finc, self.dtmax)
                    self.a *= self.fa
                self.Nsteps += 1
            else:
                self.v[:] *= 0.0
                self.a = self.astart
                self.dt *= self.fdec
                self.Nsteps = 0

        self.v += self.dt * f
        dr = self.dt * self.v
        normdr = np.sqrt(np.vdot(dr, dr))
        if normdr > self.maxmove:
            dr = self.maxmove * dr / normdr
        r = atoms.get_positions()
        atoms.set_positions(r + dr)
        self.dump((self.v, self.dt))

    def computeResidual(self):
        self.res = self.cAv - self.cvEig * self.cv
        self.dres = np.vdot(self.res, self.res) ** 0.5
        #
        # gram-schmidt for residual vector
        #
        for i in xrange(self.currentSize):
            self.dgks[i] = np.vdot(self.vlist[i, :], self.res)
            self.res -= self.dgks[i] * self.vlist[i, :]
        #
        # second gram-schmidt to make them really orthogonal
        #
        for i in xrange(self.currentSize):
            self.dgks[i] = np.vdot(self.vlist[i, :], self.res)
            self.res -= self.dgks[i] * self.vlist[i, :]
        self.resnorm = np.linalg.norm(self.res)
        self.res /= self.resnorm

        orthog = 0.0
        for i in xrange(self.currentSize):
            orthog += np.vdot(self.res, self.vlist[i]) ** 2.0
            if orthog > 1e-8:
                sys.exit("Exiting davidson procedure ... orthog = %24.16f" % orthog)
        orthog = orthog ** 0.5
        if not self.deflated:
            if VERBOSE:
                print "%3d %20.14f %20.14f  %10.4g" % (self.ciEig, self.cvEig.real, self.resnorm.real, orthog.real)
        # else:
        #    print "%3d %20.14f %20.14f %20.14f (deflated)" % (self.ciEig, self.cvEig,
        #                                                      self.resnorm, orthog)

        self.iteration += 1

    def testListFitting(self):
        true1, data1 = self.generate_data()
        true2, data2 = self.generate_data()


        model = DPMixtureModel(3, 2000, 100, 1, type='BEM')
        rs = model.fit([data1, data2])
        assert(len(rs) == 2)
        for r in rs:
            print 'mu ', r.mus
            diffs = {}
            for i in gen_mean:

                diffs[i] = np.min(np.abs(r.mus - gen_mean[i]), 0)
                # print i, gen_mean[i], diffs[i], np.vdot(diffs[i],diffs[i])
                assert(np.vdot(diffs[i], diffs[i]) < 2)

        fcm1 = FCMdata('test_fcm1', data1, ['fsc', 'ssc'], [0, 1])
        fcm2 = FCMdata('test_fcm2', data2, ['fsc', 'ssc'], [0, 1])

        c = FCMcollection('fcms', [fcm1, fcm2])

        rs = model.fit(c)
        assert(len(rs) == 2)
        for r in rs:
            print 'mu ', r.mus
            diffs = {}
            for i in gen_mean:

                diffs[i] = np.min(np.abs(r.mus - gen_mean[i]), 0)
                # print i, gen_mean[i], diffs[i], np.vdot(diffs[i],diffs[i])
                assert(np.vdot(diffs[i], diffs[i]) < 2)

def target_distrib(guess, *arg):
    """
    Keyword Arguments:
    x -- the vector of params

    *args are:
    arg[0] -- a target class object
    arg[1] -- Cl_old, fluc_lm_old
    """
    #print guess, arg
    dlm,strings,params,nl,bl = arg[0]
    Cl_old, fluc_lm_old = arg[1]
    dd = cb.update_dic(params,guess,strings)
    Cl_new = cb.generate_spectrum(dd)[:,1]
    Cl_new[:2] = 1.e-35
    #print "new = ",Cl_new[50]
    #print "old = ",Cl_old[50]
    renorm = CG.renorm_term(Cl_new,bl,nl)
    mf_lm_new = hp.almxfl(CG.generate_mfterm(dlm,Cl_new,bl,nl),renorm)
    fluc_lm_type2 = hp.almxfl(CG.generate_w1term(Cl_new,bl,nl)+CG.generate_w0term(Cl_new),renorm)
    #print "new = ",fluc_lm_type2[50]
    #print "old = ",fluc_lm_old[50]
    fluc_lm_determ = hp.almxfl(fluc_lm_old,np.sqrt(Cl_new/Cl_old))
    tt1 = -1/2.*np.real(np.vdot((dlm-hp.almxfl(mf_lm_new,bl)).T,hp.almxfl((dlm-hp.almxfl(mf_lm_new,bl)),1/nl)))
    #print tt1
    tt2 = -1/2. *np.real(np.vdot((mf_lm_new).T,hp.almxfl((mf_lm_new),1./Cl_new)))
    #print tt2
    tt3 = -1/2. *np.real(np.vdot((fluc_lm_determ).T,hp.almxfl((fluc_lm_determ),1/nl*bl**2)))
    #print tt3
    #tt4 = - 1./2 *(np.arange(1,np.size(Cl_new)+1)*np.log(Cl_new)).sum()
    ##print tt4
    return [tt1,tt2,tt3],Cl_new,fluc_lm_type2

def main():
    size=300
    A,P=setA(size,0.0+1j*0.0)
    b = np.random.rand(size,1) + 0j*np.random.rand(size,1)
    b /= np.sqrt( np.vdot(b,b) )

    xstart = np.dot(np.linalg.inv(A),b)
    xstart += 1./1.*(np.random.rand(size,1) + 0j*np.random.rand(size,1))
    condition_number = np.linalg.cond(A)
    r0=b-np.dot(A,xstart)
    res=np.vdot(r0,r0) ** 0.5
    finalx=np.dot(np.linalg.inv(A),b)
    print " ::: Making A,b matrix :::"
    print "  - condition number = %12.8f" % condition_number
    print "  - x0 residual      = %12.8f" % np.real(res)
    #for i in xrange(size):
    #    print "%20.16f %20.16f" % (np.real(finalx[i]), np.imag(finalx[i]))

    def multiplyA(vector,args):
        return np.dot(A,vector).reshape(len(vector))
    def multiplyA_precon(vector,args):
        return np.multiply(P.reshape(len(vector)),np.dot(A,vector).reshape(len(vector))).reshape(len(vector))
    gmin = gMinRes(multiplyA,b,xstart,P)

    #b = np.multiply(b,P)
    #gmin = gMinRes(multiplyA_precon,b,xstart,P)
    sol = gmin.get_solution()

    print "|Ax-b| = ", np.linalg.norm(np.dot(A,sol) - b)

def power_method(file_name, error_tol, u0):

    # Creates a matrix from the .dat file
    A = np.genfromtxt(file_name, delimiter=" ")

    m = int(A.shape[0]) #rows of A

    u = np.asarray(np.asmatrix(u0))
    uRows = int(u.shape[0]) #rows of u
    uCols = int(u.shape[1]) #columns of u

    # Sets the tolerance
    tol = error_tol

    # Sets the initial number of iterations
    iteration = 0

    # Sets an array with one entry 0 to use for the error in the while loop
    eigenvalues = [0]

    # While number of iterations is less than 100, the matrix
    # A is raised to a power equal to the iteration and multiplied
    # by the original u0 that was given as an input. The dominant
    # eigenvector is found and from that, the dominant eigenvalue
    # is found.
    while iteration < 100:
        copyA = LA.matrix_power(A, iteration+1)
        x = np.zeros(shape=(m, uCols))
        for i in range(m):
            for j in range(uCols):
                for k in range(uRows):
                    x[i][j] += (copyA[i][k] * u[k][j]) #Multiplies matrix A and u
        eigenvector = x / LA.norm(x) # Finds the dominant eigenvector
        eigenRows = int(eigenvector.shape[0]) #rows of dominant eigenvector
        eigenCols = int(eigenvector.shape[1]) #columns of eigenvector

        Ax = np.zeros(shape = (m, eigenCols))

        for i in range(m):
            for j in range(eigenCols):
                for k in range(eigenRows):
                    Ax[i][j] += A[i][k] * eigenvector[k][j]
        Axx = np.vdot(Ax, eigenvector)
        xx = np.vdot(eigenvector, eigenvector)
        eigenvalue = Axx / xx # Finds the dominant eigenvalue

        eigenvalues.append(eigenvalue)

        if (np.absolute(eigenvalues[iteration+1] - eigenvalues[iteration])) <= tol:
            break

        iteration += 1

    print "Dominant eigenvalue = ", eigenvalue
    print "Dominant eigenvector =\n", eigenvector

    if iteration >= 100:
        print "Did not converge after 100 iterations."
    else:
        print "Took " + str(iteration) + " iterations to converge."