def TestStringRepresentations():
     # Test string representations (note leading spaces)
     Julian.day_offset = mpf("0")
     Julian.interval_representation = "c"
     j1, j2 = mpf("2451545.0"), mpf("2451545.5")
     assert str(Julian(j1)) == " 1Jan2000:12:00"
     assert str(Julian(mpi(j1, j2))) == " <<1Jan2000:12:00, 2Jan2000:00:00>>"
     Julian.day_offset = mpf("0.5")
     assert str(Julian(j1)) == " 1Jan2000:00:00"
     assert str(Julian(mpi(j1, j2))) == " <<1Jan2000:00:00, 1Jan2000:12:00>>"

def LDL(mat):
    """
    Algorithm for numeric LDL factization, exploiting sparse structure.

    This function is a modification of scipy.sparse.SparseMatrix._LDL_sparse,
    allowing mpmath.mpi interval arithmetic objects as entries.

    L, D are SparseMatrix objects. However we assign values through _smat member
    to avoid type conversions to Rational.
    """
    Lrowstruc = mat.row_structure_symbolic_cholesky()
    print 'Number of entries in L: ', np.sum(map(len, Lrowstruc))
    L = SparseMatrix(mat.rows, mat.rows,
                     dict([((i, i), mpi(0)) for i in range(mat.rows)]))
    D = SparseMatrix(mat.rows, mat.cols, {})
    for i in range(len(Lrowstruc)):
        for j in Lrowstruc[i]:
            if i != j:
                L._smat[(i, j)] = mat._smat.get((i, j), mpi(0))
                summ = 0
                for p1 in Lrowstruc[i]:
                    if p1 < j:
                        for p2 in Lrowstruc[j]:
                            if p2 < j:
                                if p1 == p2:
                                    summ += L[i, p1]*L[j, p1]*D[p1, p1]
                            else:
                                break
                    else:
                        break
                L._smat[(i, j)] = L[i, j] - summ
                L._smat[(i, j)] = L[i, j] / D[j, j]

            elif i == j:
                D._smat[(i, i)] = mat._smat.get((i, i), mpi(0))
                summ = 0
                for k in Lrowstruc[i]:
                    if k < i:
                        summ += L[i, k]**2*D[k, k]
                    else:
                        break
                D._smat[(i, i)] -= summ

    return L, D

 def v(self, s):
     '''Interval numbers:  allowed forms are
         1. 'a +- b'
         2. 'a (b%)'  % sign is optional
         3. '[a, b]'
     In 1, a is the midpoint of the interval and b is the half-width.
     In 2, a is the midpoint of the interval and b is the half-width.
     In 3, the interval is indicated directly.
     '''
     e = ValueError("Improperly formed interval number '%s'" %s)
     s = s.replace(" ", "")
     if "+-" in s:
         n = [mpf(strip(i)) for i in s.split("+-")]
         return mpi(n[0] - n[1], n[0] + n[1])
     elif "(" in s:
         if s[0] == "(":  # Don't confuse with a complex number (x,y)
             return None
         if ")" not in s:
             raise e
         s = s.replace(")", "")
         percent = False
         if "%" in s:
             if s[-1] != "%":
                 raise e
             percent = True
             s = s.replace("%", "")
         a, p = [mpf(strip(i)) for i in s.split("(")]
         d = p
         if percent:
             d = a*p/mpf(100)
         return mpi(a - d, a + d)
     elif "," in s:
         if "[" not in s: raise e
         if "]" not in s: raise e
         s = s.replace("[", "")
         s = s.replace("]", "")
         n = [mpf(strip(i)) for i in s.split(",")]
         return mpi(n[0], n[1])
     else:
         return None

 def mpiTests(): # Test with mpi's
     onep2 = mpi("1", "2")
     third = Rational(1,3)
     # add
     assert no_idiff(onep2 + third, onep2 + 1/mpi(3))
     # radd
     assert no_idiff(third + onep2, onep2 + 1/mpi(3))
     # sub
     assert no_idiff(onep2 - third, onep2 - 1/mpi(3))
     # rsub
     assert no_idiff(third - onep2, -onep2 + 1/mpi(3))
     # mul
     assert no_idiff(onep2 * third, onep2 * 1/mpi(3))
     # rmul
     assert no_idiff(third * onep2, onep2 * 1/mpi(3))
     # div
     assert no_idiff(onep2 / third, onep2 /(1/mpi(3)))
     # rdiv
     assert no_idiff(third / onep2, (1/mpi(3)) / onep2)

def ParseRawData(show=False):
    '''Set show to True to have the names printed to stdout.
    '''
    def Compact(s):
        return s.replace(" ", "")
    # Locations of fields in data
    locations = {
        "name" :        (0, 55),
        "value" :       (55, 77),
        "uncertainty" : (77, 98),
    }
    s = StringIO(raw_data)
    lines = s.readlines()
    constants = {}
    fix = 1
    for line in lines:
        line = strip(line)
        if not line:
            continue
        a, b = locations["name"]
        name = strip(line[a:b])
        a, b = locations["value"]
        value = Compact(line[a:b])
        a, b = locations["uncertainty"]
        uncertainty = Compact(line[a:b])
        if fix:
            if uncertainty == "(exact)":
                uncertainty = 0
            if "..." in value:
                value = value.replace("...", "")
        try:
            x = mpf(value)
            dx = mpf(uncertainty)
            if dx == 0:
                num = x
            else:
                num = mpi(x-dx, x+dx)
            constants[name] = num
        except Exception, e:
            print name
            print "  ", str(e)

def test_special_printers():
    class IntervalPrinter(LambdaPrinter):
        """Use ``lambda`` printer but print numbers as ``mpi`` intervals. """

        def _print_Integer(self, expr):
            return "mpi('%s')" % super(IntervalPrinter, self)._print_Integer(expr)

        def _print_Rational(self, expr):
            return "mpi('%s')" % super(IntervalPrinter, self)._print_Rational(expr)

    def intervalrepr(expr):
        return IntervalPrinter().doprint(expr)

    expr = sympy.sqrt(sympy.sqrt(2) + sympy.sqrt(3)) + sympy.S(1)/2

    func0 = lambdify((), expr, modules="mpmath", printer=intervalrepr)
    func1 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter)
    func2 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter())

    mpi = type(mpmath.mpi(1, 2))

    assert isinstance(func0(), mpi)
    assert isinstance(func1(), mpi)
    assert isinstance(func2(), mpi)

def test_interval_to_mpi():
    assert Interval(0, 1).to_mpi() == mpi(0, 1)
    assert Interval(0, 1, True, False).to_mpi() == mpi(0, 1)
    assert type(Interval(0, 1).to_mpi()) == type(mpi(0, 1))

def test_interval_to_mpi():
    assert Interval(0, 1).to_mpi() == mpi(0, 1)
    assert Interval(0, 1, True, False).to_mpi() == mpi(0, 1)
    assert isinstance(Interval(0, 1).to_mpi(), type(mpi(0, 1)))

 def mpi(self):
     n = mpf(self.n)/mpf(self.d)
     return mpi(n, n)

        # Rational numbers
        Rational(3, 8) : (
            "3/8", "6/16", "0 12/32", "0-15/40", "0+18/48",
            ),
        Rational(-3, 7) : (
            "-3/7", "-6/14", "-0 12/28", "-0-15/35", "-0+18/42",
            ),
        Rational(3, -7) : (
            "-3/7", "-6/14", "-0 12/28", "-0-15/35", "-0+18/42",
            ),
    }
    # Because of a bug in mpi == and != tests, we have to test them
    # differently.
    mpi_tests = {
        # Interval numbers
        mpi(1, 3) : (
            "[1, 3]", "[1.0, 3]", "[1, 3.0]", "[1.0, 3.0]",
            "[1,3]", "[1.0,3]", "[1,3.0]", "[1.0,3.0]",
            "[   1,     3]", "[   1.0,3]", "[     1,     3.0]",
            "1.5 +- 0.5", "1.5+-0.5", "1.5+-      0.5", "1.5     +-0.5", 
            "1.5      +-      0.5", "15e-1 +- 500e-3", 
            "1.5(33.33333333333333333333%)", "1.5  (33.33333333333333333333%)",
            "1.5    (     33.33333333333333333333%)",
            "1.5    (     33.33333333333333333333     % )",
            "1.5(     33.33333333333333333333    %  )",
            "1.5(33.33333333333333333333    %  )",
            ),
    }

    n = Number()
    status = 0

    Entries assumed to be integers.
    """
    f = open(filename)
    v = np.array(eval(f.readline()), dtype=int)
    f.close()
    ind = np.array(v, dtype=int)[:, :2]
    val = np.array(v)[:, 2]
    return ss.csc_matrix((val, (ind[:, 0], ind[:, 1])))


text = open('results/ldl_mpi.txt', 'w')
for number in range(1, 21):
    positive = get_matrix('./results/CHOLMOD_permuted/permuted{}.txt'.format(number))
    size = max(positive.indices)+1
    print >>text, 'Domain {}:'.format(number)
    print >>text, 'Number of entries: ', len(positive.data)
    # change integers into mpi intervals
    print >>text, 'Entries :', np.unique(positive.data)
    positive.data = map(lambda x: mpi(x), positive.data)
    positive = positive.todok()
    sympy_positive = SparseMatrix(size, size, positive)
    L, D = LDL(sympy_positive)
    D = D._smat.values()
    delta = [x.delta/x.mid for x in D]
    print >>text, "The smallest diagonal element in LDL': ", min(D)
    print >>text, 'Ratio largest/smallest :  ', max(D)/min(D)
    print >>text, "Maximal relative delta around diagonal elements: ", max(delta)
    print >>text, '\n'
text.close()

def list_of_basis(N,weight,prec=501,tol=1e-20,sv_min=1E-1,sv_max=1E15,set_dim=None):
    r""" Returns a list of pairs (r,D) forming a basis
    """
    # First we find the smallest Discriminant for each of the components
    if(set_dim<>None and set_dim >0):
        dim=set_dim
    else:
        dim=dimension_jac_cusp_forms(int(weight+0.5),N,-1)
    basislist=dict()
    num_gotten=0
    co_tmp=dict()
    num_gotten=0
    C0=1
    RF=RealField(prec)
    if(silent>1):
        print "N=",N
        print "dim=",dim
        print "sv_min=",sv_min
        print "sv_max=",sv_max
    Aold=Matrix(RF,1)
    tol0=1E-20  #tol
    # we start with the first discriminant, then the second etc.
    Z2N=IntegerModRing(2*N)
    ZZ4N=IntegerModRing(4*N)
    for Dp in [1..max(1000,100*dim)]:
        D=-Dp # we use the dual of the Weil representation
        D4N=ZZ4N(D)
        if(not(is_square(D4N))):
            continue
        for r in my_modsqrt(D4N,N):
            # I want to make sure that P_{(D,r)} is independent from the previously computed functions
            # The only sure way to do this is to compute all submatrices (to a much smaller precision than what we want at the end)
            # The candidate is [D,r] and we need to compute the vector of [D,r,D',r']
            # for all D',r' already in the list
            ltmp1=dict()
            ltmp2=dict()
            j=0
            for [Dp,rp] in basislist.values():
                ltmp1[j]=[D,r,Dp,rp]
                ltmp2[j]=[Dp,rp,D,r]
                j=j+1
            ltmp1[j]=[D,r,D,r]
            #print "Checking: D,r,D,r=",ltmp1
            ctmp1=ps_coefficients_holomorphic_vec(N,weight,ltmp1,tol0)
            # print "ctmp1=",ctmp1
            if(j >0):
                #print "Checking: D,r,Dp,rp=",ltmp2    # Data is ok?: {0: True} 
                ctmp2=ps_coefficients_holomorphic_vec(N,weight,ltmp2,tol0)
                # print "ctmp2=",ctmp2
            #print "num_gotten=",num_gotten
            A=matrix(RF,num_gotten+1)
            # The old matrixc with the elements that are already added to the basis
            # print "Aold=\n",A,"\n"
            # print "num_gotten=",num_gotten
            # print "Aold=\n",Aold,"\n"
            for k in range(Aold.nrows()):
                for l in range(Aold.ncols()):
                    A[k,l]=Aold[k,l]
                    # endfor
                    # print "A set by old=\n",A,"\n"
                    # Add the (D',r',D,r) for each D',r' in the list
            tmp=RF(1.0)
            for l in range(num_gotten):                
                # we do not use  the scaling factor when
                # determining linear independence
                # mm=RF(abs(ltmp2[l][0]))/N4
                # tmp=RF(mm**(weight-one))
                A[num_gotten,l]=ctmp2['data'][l]*tmp
                # Add the (D,r,D',r') for each D',r' in the list
                # print "ctmp1.keys()=",ctmp1.keys()
            for l in range(num_gotten+1):
                #mm=RF(abs(ltmp1[l][2]))/4N
                #tmp=RF(mm**(weight-one))
                # print "scaled with=",tmp.n(200)
                A[l,num_gotten]=ctmp1['data'][l]*tmp
            #[d,B]=mat_inverse(A) # d=det(A) 
            #if(silent>1):
            #d=det(A)
            #print "det A = ",d
            # Now we have to determine whether we have a linearly independent set or not
            dold=mpmath.mp.dps
            mpmath.mp.dps=int(prec/3.3)
            AInt=mpmath.matrix(int(A.nrows()),int(A.ncols()))
            AMp=mpmath.matrix(int(A.nrows()),int(A.ncols()))
            if(silent>0):
                print "tol0=",tol0
            for ir in range(A.nrows()):
                for ik in range(A.ncols()):
                    AInt[ir,ik]=mpmath.mp.mpi(A[ir,ik]-tol0,A[ir,ik]+tol0)
                    AMp[ir,ik]=mpmath.mpf(A[ir,ik])

            d=mpmath.det(AMp)
            di=mpmath.mp.mpi(mpmath.mp.det(AInt))
            #for ir in range(A.nrows()):
            #    for ik in range(A.ncols()):
            #        #print "A.d=",AInt[ir,ik].delta
            if(silent>0):
                print "mpmath.mp.dps=",mpmath.mp.dps
                print "det(A)=",d
                print "det(A-as-interval)=",di
#.........这里部分代码省略.........

#.........这里部分代码省略.........
    RF=RealField(prec)
    A=matrix(RF,dim)
    kappa=weight
    fourpi=RF(4.0)*pi.n(prec)
    one=RF(1.0)
    N4=RF(4*N)
    C=dict()
    if(silent>1):
        print "v=",v
        print "dim=",dim
    lastix=0
    # First set the upper right part of A
    for j in range(dim):
        ddim=dim-j
        if(silent>1):
            print "j=",j,"ddim=",ddim," lastix=",lastix
        for k in range(0,ddim):
            # need to scale with |D|^(k+0.5)
            if(silent>1):
                print "k=",k
                print "lastix+k=",lastix+k
            mm=RF(abs(v[lastix+k][0]))/N4
            tmp=RF(mm**(weight-one))
            if(silent>1):
                print "ddim+k=",ddim+k
            A[j,j+k]=Cps[lastix+k]*tmp
            C[v[lastix+k][0],v[lastix+k][1]]=Cps[lastix+k]
        lastix=lastix+k+1
    # And add the lower triangular part to mak the matrix symmetric
    for j in range(dim):
        for k in range(0,j):
            A[j,k]=A[k,j]
    # And print the gram matrix
    res['matrix']=A
    dold=mpmath.mp.dps
    mpmath.mp.dps=int(prec/3.3)
    AInt=mpmath.matrix(int(A.nrows()),int(A.ncols()))
    AMp=mpmath.matrix(int(A.nrows()),int(A.ncols()))
    for ir in range(A.nrows()):
        for ik in range(A.ncols()):
            AInt[ir,ik]=mpmath.mpi(A[ir,ik]-tol,A[ir,ik]+tol)
            AMp[ir,ik]=mpmath.mpf(A[ir,ik])
    d=mpmath.det(AMp)
    if(silent>1):
        print "det(A-as-mpmath)=",d
    di=mpmath.det(AInt)
    if(silent>1):
        print "det(A-as-interval)=",di
    res['det']=(RF(di.a),RF(di.b))
    
    filename="PS_Gram"+stN+"-"+wt+".txt"
    if(silent>1):
        print "printing to file: "+filename
    print_matrix_to_file(A,filename,'A['+str(N)+']')
    if(silent>1):
        print "A-A.transpose()=",norm(A-A.transpose())
    B=A^-1
    #[d,B]=mat_inverse(A)
    if(silent>1):
        print "A=",A.n(100)
        print "det(A)=",di
        print "Done making inverse!"
    #res['det']=d
    res['inv']=B
    mpmath.mp.dps=dold
    filename="PS_Gram-inv"+stN+"-"+wt+".txt"        
    print_matrix_to_file(B,filename,' AI['+str(N)+']')
    # first make the filename
    s='%.1e'%tol
    filename3="PS_Coeffs"+stN+"-"+wt+"-"+s+".sobj"
    # If the file already exist we load it and append the new data
    if(silent>0):
        print "saving data to ",filename3
    try:
        f=open(filename3,"read")
    except IOError:
        if(silent>0):
            print "no file before!"
        # do nothing
    else:
        if(silent>0):
            print "file: "+filename3+" exists!"
        f.close()
        Cold=load(filename3)
        for key in Cold.keys():
            #                print"key:",key
            if(not C.has_key(key)): # then we addd it
                print"key:",key," does not exist in the new version!"
                C[key]=Cold[key]
                save(C,filename3)
    ## Save the whole thing
    filename="PS_all_gram"+stN+"-"+wt+".sobj"
    save(res,filename) 
    ## our work is comleted and we can remove the file
    try:
        os.remove(filename_work)
    except os.error:
        print "Could not remove file:",filename_work
        pass
    return res