def test_evalf():
    x = symbols('x')
    R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
    # a straight line on real axis
    r = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
    p = HolonomicFunction((1 + x)*Dx**2 + Dx, x, 0, [0, 1])
    s = '0.699525841805253'
    assert sstr(p.evalf(r)[-1]) == s
    # a traingle with vertices (0, 1+i, 2)
    r = [0.1 + 0.1*I]
    for i in range(9):
        r.append(r[-1]+0.1+0.1*I)
    for i in range(10):
        r.append(r[-1]+0.1-0.1*I)
    s = '1.07530466271334 - 0.0251200594793912*I'
    assert sstr(p.evalf(r)[-1]) == s
    p = HolonomicFunction(Dx**2 + 1, x, 0, [0, 1])
    s = '0.905546532085401 - 6.93889390390723e-18*I'
    assert sstr(p.evalf(r)[-1]) == s
    # a rectangular path (0 -> i -> 2+i -> 2)
    r = [0.1*I]
    for i in range(9):
        r.append(r[-1]+0.1*I)
    for i in range(20):
        r.append(r[-1]+0.1)
    for i in range(10):
        r.append(r[-1]-0.1*I)
    p = HolonomicFunction(Dx**2 + 1, x, 0, [1,1]).evalf(r)
    s = '0.501421652861245 - 3.88578058618805e-16*I'
    assert sstr(p[-1]) == s

def test_integrate():
    x = symbols('x')
    R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx')
    p = expr_to_holonomic(sin(x)**2/x, x0=1).integrate((x, 2, 3))
    q = '0.166270406994788'
    assert sstr(p) == q
    p = expr_to_holonomic(sin(x)).integrate((x, 0, x)).to_expr()
    q = 1 - cos(x)
    assert p == q
    p = expr_to_holonomic(sin(x)).integrate((x, 0, 3))
    q = 1 - cos(3)
    assert p == q
    p = expr_to_holonomic(sin(x)/x, x0=1).integrate((x, 1, 2))
    q = '0.659329913368450'
    assert sstr(p) == q
    p = expr_to_holonomic(sin(x)**2/x, x0=1).integrate((x, 1, 0))
    q = '-0.423690480850035'
    assert sstr(p) == q
    p = expr_to_holonomic(sin(x)/x)
    assert p.integrate(x).to_expr() == Si(x)
    assert p.integrate((x, 0, 2)) == Si(2)
    p = expr_to_holonomic(sin(x)**2/x)
    q = p.to_expr()
    assert p.integrate(x).to_expr() == q.integrate((x, 0, x))
    assert p.integrate((x, 0, 1)) == q.integrate((x, 0, 1))
    assert expr_to_holonomic(1/x).integrate(x).to_expr() == log(x)
    p = expr_to_holonomic((x + 1)**3*exp(-x), x0=-1, lenics=4).integrate(x).to_expr()
    q = (-x**3 - 6*x**2 - 15*x + 6*exp(x + 1) - 16)*exp(-x)
    assert p == q

def test_extended_domain_in_expr_to_holonomic():
    x = symbols('x')
    p = expr_to_holonomic(1.2*cos(3.1*x))
    assert p.to_expr() == 1.2*cos(3.1*x)
    assert sstr(p.integrate(x).to_expr()) == '0.387096774193548*sin(3.1*x)'
    _, Dx = DifferentialOperators(RR.old_poly_ring(x), 'Dx')
    p = expr_to_holonomic(1.1329138213*x)
    q = HolonomicFunction((-1.1329138213) + (1.1329138213*x)*Dx, x, 0, {1: [1.1329138213]})
    assert p == q
    assert p.to_expr() == 1.1329138213*x
    assert sstr(p.integrate((x, 1, 2))) == sstr((1.1329138213*x).integrate((x, 1, 2)))
    y, z = symbols('y, z')
    p = expr_to_holonomic(sin(x*y*z), x=x)
    assert p.to_expr() == sin(x*y*z)
    assert p.integrate(x).to_expr() == (-cos(x*y*z) + 1)/(y*z)
    p = expr_to_holonomic(sin(x*y + z), x=x).integrate(x).to_expr()
    q = (cos(z) - cos(x*y + z))/y
    assert p == q
    a = symbols('a')
    p = expr_to_holonomic(a*x, x)
    assert p.to_expr() == a*x
    assert p.integrate(x).to_expr() == a*x**2/2
    D_2, C_1 = symbols("D_2, C_1")
    p = expr_to_holonomic(x) + expr_to_holonomic(1.2*cos(x))
    p = p.to_expr().subs(D_2, 0)
    assert p - x - 1.2*cos(1.0*x) == 0
    p = expr_to_holonomic(x) * expr_to_holonomic(1.2*cos(x))
    p = p.to_expr().subs(C_1, 0)
    assert p - 1.2*x*cos(1.0*x) == 0

def given():
    a = 1
    b = R(randrange(3, 6))
    c = R('0.6') + randrange(5) * R(1, 10)
    e = c + randrange(4) * R(1, 10)
    d = a + randrange(5, 9)
    f = e + randrange(4, 7)
    j = f - R(1, 2)
    g = d + R(1, 2)
    h = d + randrange(4, 9)
    k = R(b - c, a * a)
    u = -k * x ** 2 + b
    v = (e - c) / (d - a) * (x - a) + c
    l = (f - e) / (h - d) ** R(1, 2)
    w = l * (x - d) ** R(1, 2) + e
    z = l * (x - d) ** R(1, 2)
    dc = Struct(
        u=sympy.sstr(u),
        v=sympy.sstr(v),
        w=sympy.sstr(w),
        z=sympy.sstr(z),
        d=d,
        g=g,
        h=h)
    return dc

    def description(self):
        from sympy import sstr

        s = 'RobotDef instance: ' + self.name + '\n'
        s += '  DH parameters - ' + str(self.dh_convention) + ' convention:\n'
        s += '    joint - alpha,  a,  d,  theta\n'
        for i, dh in enumerate(self.dh_parms):
            s += '      %2i  - ' % (i+1) + \
                sstr(dh).replace('(', '').replace(')', '') + '\n'
        s += '  gravity acceleration: ' + sstr(list(self.gravityacc)) + '^T\n'
        s += '  friction model: ' + str(self.frictionmodel) + '\n'
        return s

def test_complex_space():
    c1 = ComplexSpace(2)
    assert isinstance(c1, ComplexSpace)
    assert c1.dimension == 2
    assert sstr(c1) == "C(2)"
    assert srepr(c1) == "ComplexSpace(Integer(2))"

    n = Symbol("n")
    c2 = ComplexSpace(n)
    assert isinstance(c2, ComplexSpace)
    assert c2.dimension == n
    assert sstr(c2) == "C(n)"
    assert srepr(c2) == "ComplexSpace(Symbol('n'))"
    assert c2.subs(n, 2) == ComplexSpace(2)

    def polyrs(self, uhat):
        "Compute the polynomial for the response surface."
        import sympy

        dim = len(self.params)
        var = []
        for i, p in enumerate(self.params):
            var.append(sympy.Symbol(str(p.name)))

        chaos = chaos_sequence(dim, self.level)
        chaos = np.int_(chaos)

        for d in range(0, dim):
            poly = np.array(map(lambda x: sympy.legendre(x, var[d]), chaos[:, d]))
            if d == 0:
                s = poly
            else:
                s *= poly

        eqn = np.sum(uhat * s)
        for i, p in enumerate(self.params):
            pmin, pmax = p.pdf.srange
            c = (pmax + pmin) / 2.0
            s = (pmax - pmin) / 2.0
            eqn = eqn.subs(var[i], (var[i] - c)/s)
        return sympy.sstr(eqn.expand().evalf())

def calc(g):
    nm = 1
    for i, ae in enumerate(g.en):
        nm = nm * simplify(g.bn[i]) ** ae
    for i, ae in enumerate(g.ed):
        nm = nm / simplify(g.bd[i]) ** ae
    return [sstr(simplify(nm))]

def test_extended_domain_in_expr_to_holonomic():
    x = symbols('x')
    p = expr_to_holonomic(1.2*cos(3.1*x), domain=RR)
    assert p.to_expr() == 1.2*cos(3.1*x)
    assert sstr(p.integrate(x).to_expr()) == '0.387096774193548*sin(3.1*x)'
    _, Dx = DifferentialOperators(RR.old_poly_ring(x), 'Dx')
    p = expr_to_holonomic(1.1329138213*x, domain=RR, lenics=2)
    q = HolonomicFunction((-1.1329138213) + (1.1329138213*x)*Dx, x, 0, [0, 1.1329138213])
    assert p == q
    assert p.to_expr() == 1.1329138213*x
    assert sstr(p.integrate((x, 1, 2))) == sstr((1.1329138213*x).integrate((x, 1, 2)))
    y, z = symbols('y, z')
    p = expr_to_holonomic(sin(x*y*z), x=x, domain=ZZ[y, z])
    assert p.to_expr() == sin(x*y*z)
    assert p.integrate(x).to_expr() == (-cos(x*y*z) + 1)/(y*z)
    p = expr_to_holonomic(sin(x*y + z), x=x, domain=ZZ[y, z]).integrate(x).to_expr()
    q = (cos(z) - cos(x*y + z))/y
    assert p == q

def test_integrate():
    x = symbols('x')
    R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx')
    p = from_sympy(sin(x)**2/x, x0=1).integrate((x, 2, 3))
    q = '0.166270406994788'
    assert sstr(p) == q
    p = from_sympy(sin(x)).integrate((x, 0, x)).to_sympy()
    q = 1 - cos(x)
    assert p == q
    p = from_sympy(sin(x)).integrate((x, 0, 3))
    q = '1.98999246812687'
    assert sstr(p) == q
    p = from_sympy(sin(x)/x, x0=1).integrate((x, 1, 2))
    q = '0.659329913368450'
    assert sstr(p) == q
    p = from_sympy(sin(x)**2/x, x0=1).integrate((x, 1, 0))
    q = '-0.423690480850035'
    assert sstr(p) == q

def test_integrate():
    x = symbols('x')
    R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx')
    p = expr_to_holonomic(sin(x)**2/x, x0=1).integrate((x, 2, 3))
    q = '0.166270406994788'
    assert sstr(p) == q
    p = expr_to_holonomic(sin(x)).integrate((x, 0, x)).to_expr()
    q = 1 - cos(x)
    assert p == q
    p = expr_to_holonomic(sin(x)).integrate((x, 0, 3))
    q = 1 - cos(3)
    assert p == q
    p = expr_to_holonomic(sin(x)/x, x0=1).integrate((x, 1, 2))
    q = '0.659329913368450'
    assert sstr(p) == q
    p = expr_to_holonomic(sin(x)**2/x, x0=1).integrate((x, 1, 0))
    q = '-0.423690480850035'
    assert sstr(p) == q
    p = expr_to_holonomic(sin(x)/x)
    assert p.integrate(x).to_expr() == Si(x)
    assert p.integrate((x, 0, 2)) == Si(2)
    p = expr_to_holonomic(sin(x)**2/x)
    q = p.to_expr()
    assert p.integrate(x).to_expr() == q.integrate((x, 0, x))
    assert p.integrate((x, 0, 1)) == q.integrate((x, 0, 1))
    assert expr_to_holonomic(1/x, x0=1).integrate(x).to_expr() == log(x)
    p = expr_to_holonomic((x + 1)**3*exp(-x), x0=-1).integrate(x).to_expr()
    q = (-x**3 - 6*x**2 - 15*x + 6*exp(x + 1) - 16)*exp(-x)
    assert p == q
    p = expr_to_holonomic(cos(x)**2/x**2, y0={-2: [1, 0, -1]}).integrate(x).to_expr()
    q = -Si(2*x) - cos(x)**2/x
    assert p == q
    p = expr_to_holonomic(sqrt(x**2+x)).integrate(x).to_expr()
    q = (x**(3/2)*(2*x**2 + 3*x + 1) - x*sqrt(x + 1)*asinh(sqrt(x)))/(4*x*sqrt(x + 1))
    assert p == q
    p = expr_to_holonomic(sqrt(x**2+1)).integrate(x).to_expr()
    q = (sqrt(x**2+1)).integrate(x)
    assert (p-q).simplify() == 0
    p = expr_to_holonomic(1/x**2, y0={-2:[1, 0, 0]})
    r = expr_to_holonomic(1/x**2, lenics=3)
    assert p == r
    q = expr_to_holonomic(cos(x)**2)
    assert (r*q).integrate(x).to_expr() == -Si(2*x) - cos(x)**2/x

 def generate_equations(self):
     self.emit("EQUATION")
     self.indent()
     obs_names = self.model.observable_groups.keys()
     obs_exprs = [' + '.join('%g * s%s' % g for g in self.model.observable_groups[name]) for name in obs_names]
     for obs in zip(obs_names, obs_exprs):
         self.emit('%s = %s;' % obs)
     var_names = ['s%d' % i for i in range(len(self.model.species))]
     for (name, ode) in zip(var_names, self.model.odes):
         expression = sympy.sstr(ode)
         expression = re.sub(r'\*\*', '^', expression)
         self.emit('$%s = %s;' % (name, expression))
     self.outdent()

def test_from_hyper():
    x = symbols('x')
    R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
    p = hyper([1, 1], [S(3)/2], x**2/4)
    q = HolonomicFunction((4*x) + (5*x**2 - 8)*Dx + (x**3 - 4*x)*Dx**2, x, 1, [2*sqrt(3)*pi/9, -4*sqrt(3)*pi/27 + 4/3])
    r = from_hyper(p)
    assert r == q
    p = from_hyper(hyper([1], [S(3)/2], x**2/4))
    q = HolonomicFunction(-x + (-x**2/2 + 2)*Dx + x*Dx**2, x)
    x0 = 1
    y0 = '[sqrt(pi)*exp(1/4)*erf(1/2), -sqrt(pi)*exp(1/4)*erf(1/2)/2 + 1]'
    assert sstr(p.y0) == y0
    assert q.annihilator == p.annihilator

def norm_expr(a):
    '''
    >>> a = '1/4*e^(4x+4)'
    >>> norm_expr(a)
    'exp(4*x + 4)/4'

    '''
    try:
        from sympy.parsing.sympy_parser import (
            parse_expr,
            standard_transformations,
            convert_xor,
            implicit_multiplication_application)
        res = parse_expr(a, transformations=(
            standard_transformations + (convert_xor, implicit_multiplication_application,)))
        res = sstr(res.subs('e', E))
    except:
        res = a
    return res

def driver1(spinfile,interactionfile):
    """generates Hsave"""
    atom_list, jnums, jmats,N_atoms_uc=readfiles.readFiles(interactionfile,spinfile)
    N_atoms=len(atom_list)

    print 'N_atoms',N_atoms,'Natoms_uc',N_atoms_uc

    for atom in atom_list:
        print atom.neighbors
        print atom.interactions
    
    Hsave = calculate_dispersion(atom_list,N_atoms_uc,N_atoms,jmats,showEigs=False)   
    
    print 'driver1: complete'
    print Hsave
    
    Hamfile = open('Hsave.txt','w')
    Hamfile.write(sympy.sstr(Hsave))
    Hamfile.close()
    
    return Hsave

def NS(e, n=15, **options):
    return sstr(sympify(e).evalf(n, **options), full_prec=True)

 def __repr__(self):
     from sympy import sstr
     return '[' + ', '.join(sstr(x) for x in self.data) + ']'

def symbolic_analysis(circ, source=None, ac_enable=True, r0s=False, subs=None, outfile=None, verbose=3):
    """Attempt a symbolic solution of the circuit.

    **Parameters:**

    circ : circuit instance
        the circuit instance to be simulated.

    source : string, optional
        the ``part_id`` of the source to be used as input for the transfer
           function. If ``None``, no transfer function is evaluated.

    ac_enable : bool, optional
        take frequency dependency into consideration (default: True).

    r0s : bool, optional
        take transistors' output impedance into consideration (default: False)

    subs: dict, optional
        a dictionary of ``sympy.Symbol`` to be substituted. It makes solving the circuit
        easier. Eg. ``subs={R1:R2}`` - replace R1 with R2. It can be generated with
        :func:`parse_substitutions()`

    outfile : string, optional
        output filename - ``'stdout'`` means print to stdout, the default.

    verbose: int, optional
        verbosity level 0 (silent) to 6 (painful).

    **Returns:** 

    sol : dict
        a dictionary with the solutions.

    """
    if subs is None:
        subs = {}  # no subs by default

    if not ac_enable:
        printing.print_info_line(
            ("Starting symbolic DC analysis...", 1), verbose)
    else:
        printing.print_info_line(
            ("Starting symbolic AC analysis...", 1), verbose)

    printing.print_info_line(
        ("Building symbolic MNA, N and x...", 3), verbose, print_nl=False)
    mna, N, subs_g = generate_mna_and_N(
        circ, opts={'r0s': r0s}, ac=ac_enable, verbose=verbose)
    x = get_variables(circ)
    mna = mna[1:, 1:]
    N = N[1:, :]
    printing.print_info_line((" done.", 3), verbose)

    printing.print_info_line(
        ("Performing variable substitutions...", 5), verbose)
    mna, N = apply_substitutions(mna, N, subs)

    printing.print_info_line(("MNA matrix (reduced):", 5), verbose)
    printing.print_info_line((sympy.sstr(mna), 5), verbose)
    printing.print_info_line(("N matrix (reduced):", 5), verbose)
    printing.print_info_line((sympy.sstr(N), 5), verbose)

    printing.print_info_line(("Building equations...", 3), verbose)
    eq = []
    for i in _to_real_list(mna * x + N):
        eq.append(sympy.Eq(i, 0))

    x = _to_real_list(x)

    if verbose > 3:
        printing.print_symbolic_equations(eq)
    printing.print_info_line(("To be solved for:", 4), verbose)
    printing.print_info_line((str(x), 4), verbose)
    printing.print_info_line(("Solving...", 1), verbose)

    sol = sympy.solve(
            eq, x, manual=options.symb_sympy_manual_solver, simplify=True)

    for ks in sol.keys():
        sol.update({ks: sol[ks].subs(subs_g)})

    # sol = sol_to_dict(sol, x)

    if sol == {}:
        printing.print_warning("No solutions. Check the netlist.")
    else:
        printing.print_info_line(("Success!", 2), verbose)
        printing.print_info_line(("Results:", 1), verbose)
        if options.cli:
            printing.print_symbolic_results(sol)

    if source is not None:
        src = _symbol_factory(source.upper())
        printing.print_info_line(("Calculating small-signal symbolic transfer functions (%s))..." %
                                 (str(src),), 2), verbose, print_nl=False)
        tfs = calculate_gains(sol, src)
        printing.print_info_line(("done.", 2), verbose)
        printing.print_info_line(
            ("Small-signal symbolic transfer functions:", 1), verbose)
#.........这里部分代码省略.........

    def __repr__(self):
        from sympy import sstr

        return "[" + ", ".join(sstr(x) for x in self.data) + "]"

def symbolic_analysis(circ, source=None, ac_enable=True, r0s=False, subs=None, outfile=None, verbose=3):
    """Attempt a symbolic solution of the circuit.
    circ: the circuit instance to be simulated.
    source: the name (string) of the source to be used as input for the transfer
           function. If None, no transfer function is evaluated.
    ac_enable: take frequency dependency into consideration (default: True)
    r0s: take transistors' output impedance into consideration (default: False)
    subs: a dictionary of sympy Symbols to be substituted. It makes solving the circuit
          easier. Eg. {R1:R2} - replace R1 with R2. It can be generated with
          parse_substitutions()
    outfile: output filename ('stdout' means print to stdout).
    verbose: verbosity level 0 (silent) to 6 (painful).

    Returns: a dictionary with the solutions.
    """
    if subs is None:
        subs = {}  # no subs by default

    if not ac_enable:
        printing.print_info_line(
            ("Starting symbolic DC analysis...", 1), verbose)
    else:
        printing.print_info_line(
            ("Starting symbolic AC analysis...", 1), verbose)

    printing.print_info_line(
        ("Building symbolic MNA, N and x...", 3), verbose, print_nl=False)
    mna, N, subs_g = generate_mna_and_N(
        circ, opts={'r0s': r0s}, ac=ac_enable, verbose=verbose)
    x = get_variables(circ)
    mna = mna[1:, 1:]
    N = N[1:, :]
    printing.print_info_line((" done.", 3), verbose)

    printing.print_info_line(
        ("Performing variable substitutions...", 5), verbose)
    mna, N = apply_substitutions(mna, N, subs)

    printing.print_info_line(("MNA matrix (reduced):", 5), verbose)
    printing.print_info_line((sympy.sstr(mna), 5), verbose)
    printing.print_info_line(("N matrix (reduced):", 5), verbose)
    printing.print_info_line((sympy.sstr(N), 5), verbose)

    printing.print_info_line(("Building equations...", 3), verbose)
    eq = []
    for i in to_real_list(mna * x + N):
        eq.append(sympy.Eq(i, 0))

    x = to_real_list(x)
    if verbose > 4:
        printing.print_symbolic_equations(eq)
        print "To be solved for:"
        print x
        # print "Matrix is singular: ", (mna.det() == 0)
    # print -1.0*mna.inv()*N #too heavy
    # print sympy.solve_linear_system(mna.row_join(-N), x)
    printing.print_info_line(
        ("Performing auxiliary simplification...", 3), verbose)
    eq, x, sol_h = help_the_solver(eq, x)

    if len(eq):
        printing.print_info_line(("Symplified sytem:", 3), verbose)
        if verbose > 3:
            printing.print_symbolic_equations(eq)
        printing.print_info_line(("To be solved for:", 3), verbose)
        printing.print_info_line((str(x), 3), verbose)
        printing.print_info_line(("Solving...", 1), verbose)

        if options.symb_internal_solver:
            sol = local_solve(eq, x)
        else:
            sol = sympy.solve(
                eq, x, manual=options.symb_sympy_manual_solver, simplify=True)
        if sol is not None:
            sol.update(sol_h)
        else:
            sol = sol_h
    else:
        printing.print_info_line(
            ("Auxiliary simplification solved the problem.", 3), verbose)
        sol = sol_h

    for ks in sol.keys():
        sol.update({ks: sol[ks].subs(subs_g)})

    # sol = sol_to_dict(sol, x)

    if sol == {}:
        printing.print_warning("No solutions. Check the netlist.")
    else:
        printing.print_info_line(("Success!", 2), verbose)
        printing.print_info_line(("Results:", 1), verbose)
        if options.cli:
            printing.print_symbolic_results(sol)

    if source is not None:
        src = sympy.Symbol(source.upper(), real=True)
        printing.print_info_line(("Calculating small-signal symbolic transfer functions (%s))..." %
                                 (str(src),), 2), verbose, print_nl=False)
        tfs = calculate_gains(sol, src)
#.........这里部分代码省略.........