def main():
    print(symbols('a b c d'))
    a, b, c, d = symbols('a b c d')

    group4 = [a, b, c, d]
    print(group4)
    four = group_solve(group4)
    print(four)
    print(expand(four))
    print(four.collect(a))
    print(four.subs(a, 10))
    print()

    group4_ = add1(group4)
    print(group4_)
    four1 = group_solve(group4_)
    print(four1)
    print(expand(four1))
    print(four1.collect(a))
    print()

    print(13717421)

    for i in range(1, 10 ** 4):
        sum_power(i)

 def _solve_recur(expr, rhs=S(0)):
     if expr == var:
         return rhs
     expr = expand(expr)
     if isinstance(expr, Mul):
         lhs = S(1)
         # Try by rank
         l_coeff, var_expr, r_coeff = expr_coeff(expr, var)
         lhs = var_expr
         rhs = Inverse(l_coeff) * rhs * Inverse(r_coeff)
         return _solve_recur(lhs, rhs)
     if isinstance(expr, Add):
         lhs = 0
         for arg in expr.args:
             if var in arg:
                 lhs += arg
             else:
                 rhs -= arg
         if isinstance(lhs, Add):
             coeff = lhs.coeff(var)
             if expand(coeff * var) == lhs:
                 rhs /= coeff
                 lhs = var
         return _solve_recur(lhs, rhs)
     if isinstance(expr, Transpose):
         return _solve_recur(expr.args[0], Transpose(rhs))
     raise NotImplementedError("Can't handle expr of type %s" % type(expr))

def case0(f, N=3):
    B = 1 - x ** 3
    dBdx = sm.diff(B, x)

    # Compute basis functions and their derivatives
    phi = {0: [x ** (i + 1) * (1 - x) for i in range(N + 1)]}
    phi[1] = [sm.diff(phi_i, x) for phi_i in phi[0]]

    def integrand_lhs(phi, i, j):
        return phi[1][i] * phi[1][j]

    def integrand_rhs(phi, i):
        return f * phi[0][i] - dBdx * phi[1][i]

    Omega = [0, 1]

    u_bar = solve(integrand_lhs, integrand_rhs, phi, Omega, verbose=True, numint=False)
    u = B + u_bar
    print "solution u:", sm.simplify(sm.expand(u))

    # Calculate analytical solution

    # Solve -u''=f by integrating f twice
    f1 = sm.integrate(f, x)
    f2 = sm.integrate(f1, x)
    # Add integration constants
    C1, C2 = sm.symbols("C1 C2")
    u_e = -f2 + C1 * x + C2
    # Find C1 and C2 from the boundary conditions u(0)=0, u(1)=1
    s = sm.solve([u_e.subs(x, 0) - 1, u_e.subs(x, 1) - 0], [C1, C2])
    # Form the exact solution
    u_e = -f2 + s[C1] * x + s[C2]
    print "analytical solution:", u_e
    # print 'error:', u - u_e  # many terms - which cancel
    print "error:", sm.expand(u - u_e)

def alphaBetaToTransfer(adjMat, alphas, betas):
    s = sympy.Symbol('s')
    subs = {}

    # iterate over the columns of the adjacency matrix
    for i, col in enumerate(adjMat.getA()):
        expr = sympy.S(0)
        i += 1
        # iterate over the elements in each column of the adjacency matrix
        # this is to build the replacement expression for what were once "output to env"s
        # but are now the turnover rates of each node
        for j, elem in enumerate(col):
            j += 1
            if elem != 0:
                expr += (-sympy.var('a_%d%d' % (j,i)))

        subs[sympy.var('a_%d%d' % (i,i))] = expr

    #print "*** VAR->EXPR SUBSTITUTIONS: "
    #print subs

    for pos in range(0, len(betas)):
        orderedKeysBeta = getOrderedKeys(betas[pos])
        for key in orderedKeysBeta:
            betas[pos][s**key] = sympy.simplify(sympy.expand(betas[pos][s**key].xreplace(subs)))

        #finished all the betas so add the alphas, but only once
        orderedKeysAlpha = getOrderedKeys(alphas[pos])
        for key in orderedKeysAlpha:
            alphas[pos][s**key] = sympy.simplify(sympy.expand(alphas[pos][s**key].xreplace(subs)))

def test_one_dof():
    # This is for a 1 dof spring-mass-damper case.
    # It is described in more detail in the Kane docstring.
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, c, k = symbols('m c k')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, u * N.x)

    kd = [qd - u]
    FL = [(P, (-k * q - c * u) * N.x)]
    pa = Particle()
    pa.mass = m
    pa.point = P
    BL = [pa]

    KM = Kane(N)
    KM.coords([q])
    KM.speeds([u])
    KM.kindiffeq(kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    assert expand(rhs[0]) == expand(-(q * k + u * c) / m)
    assert KM.linearize() == (Matrix([[0, 1], [k, c]]), Matrix([]))

def _interpolate_Function(expr):
    path = 'None'
    factor = 'None'
    change = False
    summand_0 = 'None'
    summand_1 = 'None'
    res = expr
    for i in np.arange(len(expr.args)):
        argument = sympy.expand(expr.args[i])
        if argument.func == sympy.Add:
            for j in np.arange(len(argument.args)):
                summand = argument.args[j]
                if summand.func == sympy.Mul:
                    for k in np.arange(len(summand.args)):
                        temp = 0
                        if summand.args[k] == sympy.Symbol('a'):
                            temp = sympy.Mul(sympy.Mul(*summand.args[:k]),
                                             sympy.Mul(*summand.args[k+1:]))
                            #print(temp)
                        if not temp == int(temp):
                            #print('found one')
                            factor = (temp)
                            path = np.array([i, j, k])
    if not factor == 'None':
        change = True

        sign = np.sign(factor)
        offsets = np.array([int(factor), sign * (int(sign * factor) + 1)])
        weights = 1/np.abs(offsets - factor)
        weights = weights/np.sum(weights)

        res = (  weights[0] * _interchange(expr, offsets[0] * sympy.Symbol('a'), path[:-1])
               + weights[1] * _interchange(expr, offsets[1] * sympy.Symbol('a'), path[:-1]))

    return sympy.expand(res), change

def _interpolate_expression(expr):
    change = False
    expr = sympy.expand(expr)
    res = expr

    if isinstance(expr, sympy.Function):# and not change:
        path = np.array([])
        temp, change = _interpolate_Function(_follow_path(expr, path))
        res = _interchange(expr, temp, path)

    for i in np.arange(len(expr.args)):
        path = np.array([i])
        if isinstance(_follow_path(expr, path), sympy.Function) and not change:
            temp, change = _interpolate_Function(_follow_path(expr, path))
            res = _interchange(expr, temp, path)

        for j in np.arange(len(expr.args[i].args)):
            path = np.array([i, j])
            if isinstance(_follow_path(expr, path), sympy.Function) and not change:
                temp, change = _interpolate_Function(_follow_path(expr, path))
                res = _interchange(expr, temp, path)

    if change:
        res = _interpolate_expression(res)

    return sympy.expand(res)

def test_SVCVS_laplace_d3_n1():
    """Test VCCS with a laplace defined transfer function with second order
    numerator and third order denominator
    """

    pycircuit.circuit.circuit.default_toolkit = symbolic
    cir = SubCircuit()

    n1,n2 = cir.add_nodes('1','2')

    b0,a0,a1,a2,a3,Gdc = [sympy.Symbol(symname, real=True) for
                                symname in 'b0,a0,a1,a2,a3,Gdc'
                                .split(',')]

    s = sympy.Symbol('s', complex=True)

    cir['VS']   = VS( n1, gnd, vac=1)
    cir['VCVS'] = SVCVS( n1, gnd, n2, gnd,
                        denominator = [a0, a1, a2, a3],
                        numerator   = [b0, 0, 0])

    res = AC(cir, toolkit=symbolic).solve(s, complexfreq=True)

    assert_equal(sympy.cancel(sympy.expand(res.v(n2,gnd))),
                 sympy.expand((b0*s*s)/(a0*s*s*s+a1*s*s+a2*s+a3)))

def tests():
    x, y, z = symbols('x,y,z')
    #print(x + x + 1)
    expr = x**2 - y**2
    factors = factor(expr)
    print(factors, " | ", expand(factors))
    pprint(expand(factors))

def test_one_dof():
    # This is for a 1 dof spring-mass-damper case.
    # It is described in more detail in the KanesMethod docstring.
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, c, k = symbols('m c k')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, u * N.x)

    kd = [qd - u]
    FL = [(P, (-k * q - c * u) * N.x)]
    pa = Particle('pa', P, m)
    BL = [pa]

    KM = KanesMethod(N, [q], [u], kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    assert expand(rhs[0]) == expand(-(q * k + u * c) / m)
    assert (KM.linearize(A_and_B=True, new_method=True)[0] ==
            Matrix([[0, 1], [-k/m, -c/m]]))

    # Ensure that the old linearizer still works and that the new linearizer
    # gives the same results. The old linearizer is deprecated and should be
    # removed in >= 0.7.7.
    M_old = KM.mass_matrix_full
    # The old linearizer raises a deprecation warning, so catch it here so
    # it doesn't cause py.test to fail.
    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=SymPyDeprecationWarning)
        F_A_old, F_B_old, r_old = KM.linearize()
    M_new, F_A_new, F_B_new, r_new = KM.linearize(new_method=True)
    assert simplify(M_new.inv() * F_A_new - M_old.inv() * F_A_old) == zeros(2)

def test_pend():
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, l, g = symbols('m l g')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, -l * u * sin(q) * N.x + l * u * cos(q) * N.y)
    kd = [qd - u]

    FL = [(P, m * g * N.x)]
    pa = Particle()
    pa.mass = m
    pa.point = P
    BL = [pa]

    KM = Kane(N)
    KM.coords([q])
    KM.speeds([u])
    KM.kindiffeq(kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    rhs.simplify()
    assert expand(rhs[0]) == expand(-g / l * sin(q))

def T_exact_symbolic(verbose=False):
    """Compute the exact solution formula via sympy."""
    # sol1: solution for t < t_star,
    # sol2: solution for t > t_star
    import sympy as sym
    T0 = sym.symbols('T0')
    k = sym.symbols('k', positive=True)
    # Piecewise linear T_sunction
    t, t_star, C0, C1 = sym.symbols('t t_star C0 C1')
    T_s = C0
    I = sym.integrate(sym.exp(k*t)*T_s, (t, 0, t))
    sol1 = T0*sym.exp(-k*t) + k*sym.exp(-k*t)*I
    sol1 = sym.simplify(sym.expand(sol1))
    if verbose:
        # Some debugging print
        print 'solution t < t_star:', sol1
        #print sym.latex(sol1)
    T_s = C1
    I = sym.integrate(sym.exp(k*t)*C0, (t, 0, t_star)) + \
        sym.integrate(sym.exp(k*t)*C1, (t, t_star, t))
    sol2 = T0*sym.exp(-k*t) + k*sym.exp(-k*t)*I
    sol2 = sym.simplify(sym.expand(sol2))
    if verbose:
        print 'solution t > t_star:', sol2
        #print sym.latex(sol2)

    # Convert to numerical functions
    exact0 = sym.lambdify([t, C0, k, T0],
                          sol1, modules='numpy')
    exact1 = sym.lambdify([t, C0, C1, t_star, k, T0],
                          sol2, modules='numpy')
    return exact0, exact1

 def gen_errs(self, a, b, n, x, y):
     """
     What is the best thing to do here?
     """
     errs = set()
     
     # Generate som obvious candidates
     errs.add(sym.expand(self.gen_prob(-a, -b, n, x, y)))
     if a < 0:
         errs.add(sym.expand(self.gen_prob(-a, b, n, x, y)))
     if b < 0:
         errs.add(sym.expand(self.gen_prob(a, -b, n, x, y)))
    
     
     expr = sym.expand(self.gen_prob(a, b, n, x, y))
     coeffs = sym.Poly(expr).coeffs()
     
     for i in range(4):
         
         expr1 = self.gen_poly(map(lambda l: random.choice([-1, 1]) * l, coeffs), x, y)
         expr2 = self.gen_poly(map(lambda l: random.choice([-2, -1, 0, 1, 2]) * l, coeffs), x, y)
         errs.update([expr1, expr2, expr - (expr2 - sym.LM(expr))])
     
     errs = list(errs)
     errs_ = [err for err in errs if err != expr]
     random.shuffle(errs_)
     return errs

def test_two_dof():
    # This is for a 2 d.o.f., 2 particle spring-mass-damper.
    # The first coordinate is the displacement of the first particle, and the
    # second is the relative displacement between the first and second
    # particles. Speeds are defined as the time derivatives of the particles.
    q1, q2, u1, u2 = dynamicsymbols('q1 q2 u1 u2')
    q1d, q2d, u1d, u2d = dynamicsymbols('q1 q2 u1 u2', 1)
    m, c1, c2, k1, k2 = symbols('m c1 c2 k1 k2')
    N = ReferenceFrame('N')
    P1 = Point('P1')
    P2 = Point('P2')
    P1.set_vel(N, u1 * N.x)
    P2.set_vel(N, (u1 + u2) * N.x)
    kd = [q1d - u1, q2d - u2]

    # Now we create the list of forces, then assign properties to each
    # particle, then create a list of all particles.
    FL = [(P1, (-k1 * q1 - c1 * u1 + k2 * q2 + c2 * u2) * N.x), (P2, (-k2 *
        q2 - c2 * u2) * N.x)]
    pa1 = Particle('pa1', P1, m)
    pa2 = Particle('pa2', P2, m)
    BL = [pa1, pa2]

    # Finally we create the KanesMethod object, specify the inertial frame,
    # pass relevant information, and form Fr & Fr*. Then we calculate the mass
    # matrix and forcing terms, and finally solve for the udots.
    KM = KanesMethod(N, q_ind=[q1, q2], u_ind=[u1, u2], kd_eqs=kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    assert expand(rhs[0]) == expand((-k1 * q1 - c1 * u1 + k2 * q2 + c2 * u2)/m)
    assert expand(rhs[1]) == expand((k1 * q1 + c1 * u1 - 2 * k2 * q2 - 2 *
                                    c2 * u2) / m)

    def test_apply_substitutions(self):

        def apply_correct_substitutions(monomial, substitutions):
            if isinstance(monomial, int) or isinstance(monomial, float):
                return monomial
            original_monomial = monomial
            changed = True
            while changed:
                for lhs, rhs in substitutions.items():
                    monomial = monomial.subs(lhs, rhs)
                if original_monomial == monomial:
                    changed = False
                original_monomial = monomial
            return monomial

        length, h, U, t = 2, 3.8, -6, 1
        fu = generate_operators('fu', length)
        fd = generate_operators('fd', length)
        _b = flatten([fu, fd])
        hamiltonian = 0
        for j in range(length):
            hamiltonian += U * (Dagger(fu[j])*Dagger(fd[j]) * fd[j]*fu[j])
            hamiltonian += -h/2*(Dagger(fu[j])*fu[j] - Dagger(fd[j])*fd[j])
            for k in get_neighbors(j, len(fu), width=1):
                hamiltonian += -t*Dagger(fu[j])*fu[k]-t*Dagger(fu[k])*fu[j]
                hamiltonian += -t*Dagger(fd[j])*fd[k]-t*Dagger(fd[k])*fd[j]
        substitutions = fermionic_constraints(_b)
        monomials = expand(hamiltonian).as_coeff_mul()[1][0].as_coeff_add()[1]
        substituted_hamiltonian = sum([apply_substitutions(monomial,
                                                           substitutions)
                                       for monomial in monomials])
        correct_hamiltonian = sum([apply_correct_substitutions(monomial,
                                                               substitutions)
                                   for monomial in monomials])
        self.assertTrue(substituted_hamiltonian == expand(correct_hamiltonian))

 def test_fast_substitute(self):
     f = generate_variables(2, name="f")
     substitutions = {}
     substitutions[Dagger(f[0]) * f[0]] = -f[0] * Dagger(f[0])
     monomial = Dagger(f[0]) * f[0]
     lhs = Dagger(f[0]) * f[0]
     rhs = -f[0] * Dagger(f[0])
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial.subs(lhs, rhs))
     monomial = Dagger(f[0]) * f[0] ** 2
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial.subs(lhs, rhs))
     monomial = Dagger(f[0]) ** 2 * f[0]
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial.subs(lhs, rhs))
     monomial = Dagger(f[0]) ** 2 * f[0]
     lhs = Dagger(f[0]) ** 2
     rhs = -f[0] * Dagger(f[0])
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial.subs(lhs, rhs))
     g = generate_variables(2, name="g")
     monomial = 2 * g[0] ** 3 * g[1] * Dagger(f[0]) ** 2 * f[0]
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial.subs(lhs, rhs))
     monomial = S.One
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial.subs(lhs, rhs))
     monomial = 5
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial)
     monomial = 2 * g[0] ** 3 * g[1] * Dagger(f[0]) ** 2 * f[0] + f[1]
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == monomial.subs(lhs, rhs))
     monomial = f[1] * Dagger(f[0]) ** 2 * f[0]
     lhs = f[1]
     rhs = 1.0 + f[0]
     self.assertTrue(fast_substitute(monomial, lhs, rhs) == expand(monomial.subs(lhs, rhs)))
     monomial = f[1] ** 2 * Dagger(f[0]) ** 2 * f[0]
     result = fast_substitute(fast_substitute(monomial, lhs, rhs), lhs, rhs)
     self.assertTrue(result == expand(monomial.subs(lhs, rhs)))

def test_order_noncommutative():
    A = Symbol("A", commutative=False)
    assert Order(A + A * x, x) == Order(1, x)
    assert (A + A * x) * Order(x) == Order(x)
    assert (A * x) * Order(x) == Order(x ** 2, x)
    assert expand((1 + Order(x)) * A * A * x) == A * A * x + Order(x ** 2, x)
    assert expand((A * A + Order(x)) * x) == A * A * x + Order(x ** 2, x)
    assert expand((A + Order(x)) * A * x) == A * A * x + Order(x ** 2, x)

def test_piecewise_fold_expand():
    p1 = Piecewise((1,Interval(0,1,False,True)),(0,True))

    p2 = piecewise_fold(expand((1-x)*p1))
    assert p2 == Piecewise((1 - x, Interval(0,1,False,True)), \
        (Piecewise((-x, Interval(0,1,False,True)), (0, True)), True))

    p2 = expand(piecewise_fold((1-x)*p1))
    assert p2 == Piecewise((1 - x, Interval(0,1,False,True)), (0, True))

def test_meijerint():
    from sympy import symbols, expand, arg

    s, t, mu = symbols("s t mu", real=True)
    assert integrate(
        meijerg([], [], [0], [], s * t) * meijerg([], [], [mu / 2], [-mu / 2], t ** 2 / 4), (t, 0, oo)
    ).is_Piecewise
    s = symbols("s", positive=True)
    assert integrate(x ** s * meijerg([[], []], [[0], []], x), (x, 0, oo)) == gamma(s + 1)
    assert integrate(x ** s * meijerg([[], []], [[0], []], x), (x, 0, oo), meijerg=True) == gamma(s + 1)
    assert isinstance(integrate(x ** s * meijerg([[], []], [[0], []], x), (x, 0, oo), meijerg=False), Integral)

    assert meijerint_indefinite(exp(x), x) == exp(x)

    # TODO what simplifications should be done automatically?
    # This tests "extra case" for antecedents_1.
    a, b = symbols("a b", positive=True)
    assert simplify(meijerint_definite(x ** a, x, 0, b)[0]) == b ** (a + 1) / (a + 1)

    # This tests various conditions and expansions:
    meijerint_definite((x + 1) ** 3 * exp(-x), x, 0, oo) == (16, True)

    # Again, how about simplifications?
    sigma, mu = symbols("sigma mu", positive=True)
    i, c = meijerint_definite(exp(-((x - mu) / (2 * sigma)) ** 2), x, 0, oo)
    assert simplify(i) == sqrt(pi) * sigma * (erf(mu / (2 * sigma)) + 1)
    assert c is True

    i, _ = meijerint_definite(exp(-mu * x) * exp(sigma * x), x, 0, oo)
    # TODO it would be nice to test the condition
    assert simplify(i) == 1 / (mu - sigma)

    # Test substitutions to change limits
    assert meijerint_definite(exp(x), x, -oo, 2) == (exp(2), True)
    assert expand(meijerint_definite(exp(x), x, 0, I)[0]) == exp(I) - 1
    assert expand(meijerint_definite(exp(-x), x, 0, x)[0]) == 1 - exp(-exp(I * arg(x)) * abs(x))

    # Test -oo to oo
    assert meijerint_definite(exp(-x ** 2), x, -oo, oo) == (sqrt(pi), True)
    assert meijerint_definite(exp(-abs(x)), x, -oo, oo) == (2, True)
    assert meijerint_definite(exp(-(2 * x - 3) ** 2), x, -oo, oo) == (sqrt(pi) / 2, True)
    assert meijerint_definite(exp(-abs(2 * x - 3)), x, -oo, oo) == (1, True)
    assert meijerint_definite(exp(-((x - mu) / sigma) ** 2 / 2) / sqrt(2 * pi * sigma ** 2), x, -oo, oo) == (1, True)

    # Test one of the extra conditions for 2 g-functinos
    assert meijerint_definite(exp(-x) * sin(x), x, 0, oo) == (S(1) / 2, True)

    # Test a bug
    def res(n):
        return (1 / (1 + x ** 2)).diff(x, n).subs(x, 1) * (-1) ** n

    for n in range(6):
        assert integrate(exp(-x) * sin(x) * x ** n, (x, 0, oo), meijerg=True) == res(n)

    # Test trigexpand:
    assert integrate(exp(-x) * sin(x + a), (x, 0, oo), meijerg=True) == sin(a) / 2 + cos(a) / 2

 def _w_diff_dcm(self, otherframe):
     """Angular velocity from time differentiating the DCM. """
     from sympy.physics.vector.functions import dynamicsymbols
     dcm2diff = self.dcm(otherframe)
     diffed = dcm2diff.diff(dynamicsymbols._t)
     angvelmat = diffed * dcm2diff.T
     w1 = trigsimp(expand(angvelmat[7]), recursive=True)
     w2 = trigsimp(expand(angvelmat[2]), recursive=True)
     w3 = trigsimp(expand(angvelmat[3]), recursive=True)
     return -Vector([(Matrix([w1, w2, w3]), self)])