def main():

    print 'Initial metric:'
    pprint(gdd)
    print '-'*40
    print 'Christoffel symbols:'
    for i in [0,1,2,3]:
    	for k in [0,1,2,3]:
    		for l in [0,1,2,3]:
    			if Gamma.udd(i,k,l) != 0 :
		    		pprint_Gamma_udd(i,k,l)
    print'-'*40
    print'Ricci tensor:'
    for i in [0,1,2,3]:
    	for j in [0,1,2,3]:
    		if Rmn.dd(i,j) !=0:
    			pprint_Rmn_dd(i,j)
    
    print '-'*40
    #Solving EFE for A and B
    s = ( Rmn.dd(1,1)/ A(r) ) + ( Rmn.dd(0,0)/ B(r) )
    pprint (s)
    t = dsolve(s, A(r))
    pprint(t)
    metric = gdd.subs(A(r), t)
    print "metric:"
    pprint(metric)
    r22 = Rmn.dd(3,3).subs( A(r), 1/B(r))
    h = dsolve( r22, B(r) )
    pprint(h)

def test_ode10():
    f = Function("f")
    #type:2nd order, constant coefficients (two real different roots)
    eq = Eq(f(x).diff(x,x) - 3*diff(f(x),x) + 2*f(x), 0)
    assert dsolve(eq, f(x)) in [
        Symbol("C1")*exp(2*x) + Symbol("C2")*exp(x),
        Symbol("C1")*exp(x) + Symbol("C2")*exp(2*x),
    ]
    checksol(eq, f(x), dsolve(eq, f(x)))

def edo_main():
    # const = [0]*6 #a0, a1, a2, a3, a4, a5   ; nessa ordem
    print "Valores", const

    ##Criando uma variavel para guardar valores digitados anteriormente,este só vai calcular as coisas novamente se
    # e somente se parâmetros de entrada diferentes forem colocados
    constAnterior = const

    # a5, a4, a3, a2, a1, a0, xT = input_coefs()
    ###Dif equation solver
    ##Sets if it is of homogenous or inhomogenous type and type of resolution method
    if Respostas[7] == 0:
        solvedEq = dsolve(sympify(Respostas[8]), y(t), hint="nth_linear_constant_coeff_homogeneous")
        # elif (a3 != 0) or (a4 != 0) or (a5 != 0):
        # solvedEq = dsolve(sympify(eq),y(t),hint='nth_linear_constant_coeff_variation_of_parameters')

    else:
        solvedEq = dsolve(sympify(Respostas[8]), y(t), hint="nth_linear_constant_coeff_undetermined_coefficients")

    ##Transformação do tipo sympy_unity para o sympy_mul (mais operações permitidas)
    sepEq = solvedEq._args[1]
    sepEq = sepEq.evalf(prec)

    if const[5] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0), (C3, 0), (C4, 0), (C5, 0)])
    elif const[4] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0), (C3, 0), (C4, 0)])
    elif const[3] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0), (C3, 0)])
    elif const[2] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0)])
    elif const[1] != 0:
        RespPart = sepEq.subs(C1, 0)

    ##Resposta transitória alocada em  RespTran, natural em RespNat

    Respostas[3] = RespPart.evalf(prec)
    formaNatural = sepEq.subs(RespPart, 0)
    Respostas[1] = formaNatural  # Adicionando Forma natural de resposta na lista de respostas
    ## fN é a mesma coisa, mas usado por um bug bizarro do Sympy que exige uma variável sem alocações prévias quando diferenciando
    ##isso é válido no método conds_iniciais_aplicadas
    fN = formaNatural
    rP = RespPart.evalf(prec)
    raizes()
    conds_iniciais_aplicadas(fN, rP)

    respForc = Respostas[4] + Respostas[3]  # Yf = Yt + Yp
    Respostas[5] = respForc.evalf(prec)  # Adiciona Resposta Forcada a lista de respostas

    respComp = Respostas[2]  # Resposta completa p/ eqs. homogeneas
    if Respostas[7] != 0:  # Eqs. nao homogeneas
        respComp = Respostas[2] + Respostas[5]  # Respsota completa p/ eqs. nao-homogeneas

    Respostas[6] = respComp.evalf(prec)  # Adiciona Resposta Completa a lista de respostas

    print "Respostas em ODE_main", Respostas

def main():
    x = Symbol("x")
    f = Function("f")

    eq = Eq(f(x).diff(x), f(x))
    print "Solution for ", eq, " : ", dsolve(eq, f(x))

    eq = Eq(f(x).diff(x, 2), -f(x))
    print "Solution for ", eq, " : ", dsolve(eq, f(x))

    eq = Eq(x ** 2 * f(x).diff(x), -3 * x * f(x) + sin(x) / x)
    print "Solution for ", eq, " : ", dsolve(eq, f(x))

def symbolic_solve():
    """
    Solves the problem symbolically using sympy
    """
    f = Function('f')
    sol = dsolve(2 * Derivative(f(x), x, x) - 1, f(x))
    print sol

def separable_equation(g, h, hf = None):
    """
    dy/dx = g(x)*h(y)
    dy/h(y) = g(x)*dx
    """
    dy, dx = symbols('dy, dx')
    print '\nODE to solve:'
    pprint(Eq(dy/dx, g*h))
    pprint(Eq(dy/h, g*dx))
    print '\nintegrate both sides:'
    LHS, RHS = symbols('LHS, RHS')
    pprint(Eq(LHS, Integral(1/h,y)))
    H = integrate(1/h,y)
    pprint(Eq(LHS, H))
    pprint(Eq(RHS,Integral(g,x)))
    G = integrate(g, x)
    pprint(Eq(RHS, G))
    C = symbols('C')
    print '\nsolving LHS = RHS + C...'
    eq = Eq(H,G+C)
    pprint(eq)
    pprint(solve(eq, y))
    if hf:
        print '\nsolving ODE directly ...'
        pprint(dsolve(f(x).diff(x)-g*hf, f(x)))

    def __call__(self, equations, variables=None):
        if variables is None:
            variables = {}

        if equations.is_stochastic:
            raise ValueError('Cannot solve stochastic equations with this state updater')

        diff_eqs = equations.substituted_expressions

        t = Symbol('t', real=True, positive=True)
        dt = Symbol('dt', real=True, positive=True)
        t0 = Symbol('t0', real=True, positive=True)
        f0 = Symbol('f0', real=True)
        # TODO: Shortcut for simple linear equations? Is all this effort really
        #       worth it?

        code = []
        for name, expression in diff_eqs:
            rhs = expression.sympy_expr
            non_constant = _non_constant_symbols(rhs.atoms(),
                                                 variables) - set([name])
            if len(non_constant):
                raise ValueError(('Equation for %s referred to non-constant '
                                  'variables %s') % (name, str(non_constant)))
            # We have to be careful and use the real=True assumption as well,
            # otherwise sympy doesn't consider the symbol a match to the content
            # of the equation
            var = Symbol(name, real=True)
            f = sp.Function(name)
            rhs = rhs.subs(var, f(t))
            derivative = sp.Derivative(f(t), t)
            diff_eq = sp.Eq(derivative, rhs)
            general_solution = sp.dsolve(diff_eq, f(t))
            # Check whether this is an explicit solution
            if not getattr(general_solution, 'lhs', None) == f(t):
                raise ValueError('Cannot explicitly solve: ' + str(diff_eq))
            # Solve for C1 (assuming "var" as the initial value and "t0" as time)
            if Symbol('C1') in general_solution:
                if Symbol('C2') in general_solution:
                    raise ValueError('Too many constants in solution: %s' % str(general_solution))
                constant_solution = sp.solve(general_solution, Symbol('C1'))
                if len(constant_solution) != 1:
                    raise ValueError(("Couldn't solve for the constant "
                                      "C1 in : %s ") % str(general_solution))
                constant = constant_solution[0].subs(t, t0).subs(f(t0), var)
                solution = general_solution.rhs.subs('C1', constant)
            else:
                solution = general_solution.rhs.subs(t, t0).subs(f(t0), var)
            # Evaluate the expression for one timestep
            solution = solution.subs(t, t + dt).subs(t0, t)
            # only try symplifying it -- it sometimes raises an error
            try:
                solution = solution.simplify()
            except ValueError:
                pass

            code.append(name + ' = ' + sympy_to_str(solution))

        return '\n'.join(code)

def non_homo_linear(rhs, *cds):
    char_func, eq, rs = homo_linear(*cds)
    print('\nnon-homogeneous ODE:')
    pprint(Eq(eq, rhs))
    eq -= rhs
    print('\nsolving non-homogeneous linear equation...\nresult:')
    rs = dsolve(eq)
    pprint(rs)
    return char_func, eq, rs

def test_make_ode_01():
    ode, params = _make_ode_01()
    t, y, y0, k = params
    result = dsolve(ode, y[1](t))
    eq_assumption = sympy.Q.is_true(Eq(k[1], k[0]))
    refined = result.refine(~eq_assumption)
    ignore = k + y0 + (t,)
    int_const = [fs for fs in refined.free_symbols if fs not in ignore][0]
    ref = int_const*exp(-k[1]*t) - exp(-k[0]*t)*k[0]*y0[0]/(k[0] - k[1])
    assert (refined.rhs - ref).simplify() == 0

    def solve_model(self, ics, d_ics):
        model_soln = [dsolve(eq) for eq in self.model]

        def solve_constants(eq, ics, d_ics):
            udiff = Eq(d_ics[0][1], eq.rhs.diff(t))
            system = [eq.subs(ics), udiff.subs(t, 0)]
            consts = solve(system, [C1, C2])
            return eq.subs(consts)

        model_soln_f = [solve_constants(eq[1], ics[eq[0]], d_ics[eq[0]])
                        for eq in enumerate(model_soln[:len(ics)])]

        self.model_response = model_soln_f
        self.x = lambdify(t, model_soln_f[0].rhs, 'numpy')
        self.y = lambdify(t, model_soln_f[1].rhs, 'numpy')

def codegen(expr, lang, indent='    ', ics=None):
    if lang == 'C':
        code = ccode
    elif lang == 'Fortran':
        code = fcode
    else:
        raise ValueError("Lang must be 'C' or 'Fortran'")

    try:
        sol = dsolve(expr, ics=ics)
    except ValueError:
        # Not an ODE
        return code(expr)

    return ccode(sol.rhs, assign_to=sol.lhs.func.__name__)

def solveSystem(L,coordinates,initialConditions):
    """
    This is currently not in use; I'm not sure it ever will be used.    
    
    L is a sympy expression for the Lagrangian.
    coordinates is a list of tuples, like so:
    [(x,xdot),(y,ydot)]
    initialConditions are a list of the above form.
    """    
    (eulerLagrange,coordinates_t,t)=calcEL(L, coordinates)
    for (i,(EL_i,coordinates_t_i)) in\
            enumerate(zip(eulerLagrange,coordinates_t)):
        eqn=sp.dsolve(EL_i,coordinates_t_i(t))
        freeVars=filter(lambda x:x!=t,eqn.atoms())
        newFreeVars=[sp.Symbol(str(freeVar)+"_%i"%i)
            for freeVar in freeVars]

def homo_linear(*cds):
    char_func = 0
    eq = 0
    for cd in cds:
        try:
            c, d = cd
        except TypeError:
            c, d = 1, cd
        char_func += c * y**d
        eq += c * f(x).diff(x, d)
    print('\nhomogeneous ODE:')
    pprint(Eq(eq, 0))
    print('\nhomogeneous characteristic function:')
    pprint(char_func)
    print('\nsolving characteristic function...\nresult:')
    pprint(solve(char_func))
    print('\nsolving homogeneous linear equation...\nresult:')
    rs = dsolve(eq)
    pprint(rs)
    return char_func, eq, rs

def main():

    print("Initial metric:")
    pprint(gdd)
    print("-"*40)
    print("Christoffel symbols:")
    pprint_Gamma_udd(0, 1, 0)
    pprint_Gamma_udd(0, 0, 1)
    print()
    pprint_Gamma_udd(1, 0, 0)
    pprint_Gamma_udd(1, 1, 1)
    pprint_Gamma_udd(1, 2, 2)
    pprint_Gamma_udd(1, 3, 3)
    print()
    pprint_Gamma_udd(2, 2, 1)
    pprint_Gamma_udd(2, 1, 2)
    pprint_Gamma_udd(2, 3, 3)
    print()
    pprint_Gamma_udd(3, 2, 3)
    pprint_Gamma_udd(3, 3, 2)
    pprint_Gamma_udd(3, 1, 3)
    pprint_Gamma_udd(3, 3, 1)
    print("-"*40)
    print("Ricci tensor:")
    pprint_Rmn_dd(0, 0)
    e = Rmn.dd(1, 1)
    pprint_Rmn_dd(1, 1)
    pprint_Rmn_dd(2, 2)
    pprint_Rmn_dd(3, 3)
    # print()
    # print "scalar curvature:"
    # print curvature(Rmn)
    print("-"*40)
    print("Solve Einstein's equations:")
    e = e.subs(nu(r), -lam(r)).doit()
    l = dsolve(e, lam(r))
    pprint(l)
    lamsol = solve(l, lam(r))[0]
    metric = gdd.subs(lam(r), lamsol).subs(nu(r), -lamsol)  # .combine()
    print("metric:")
    pprint(metric)

 def process(self):
     """
     Procesamos la ecuacion y obtenemos los resultados
     """
     
     # Guardamos expresion
     expr = ""
     
     # Obtenemos los symbolos
     x = sp.Symbol("x")
     y = sp.Function("y")
     
     # Obtenemos la expresion
     ec = sp.sympify(self.ec)
     
     # Valor inicial (PVI)
     pvi = { y(self.x): self.y }
     
     # preparamos la EDO
     edo = sp.Eq(y(x).diff(x), ec)
     expr += "EDO:\n\t"+str(edo)
     
     # Despejamos Y
     res = sp.dsolve(y(x).diff(x) - ec)
     # Obtenemos y(x) = f(x)
     expr += "\nEDO resuelta:\n\t"+str(res)
     
     # reemplazamos PVI
     c_eq = sp.Eq(res.lhs.subs(x, 0).subs(pvi), res.rhs.subs(x, 0))
     expr += "\nRemplazamos PVI:\n\t"+str(c_eq)
     
     # Obtenemos el valor de la constante
     c = sp.solve(c_eq)
     expr += "\nValor de C1:\n\t"+str(c[0])
     
     # almacenamos los valores de interes
     self.c = c[0] # valor de C1
     self.res = res # ecuacion resuelta
     
     # retornamos el resultado
     return expr

def main():

    print "Initial metric:"
    pprint(gdd)
    print "-"*40
    print "Christoffel symbols:"
    pprint_Gamma_udd(0,1,0)
    pprint_Gamma_udd(0,0,1)
    print
    pprint_Gamma_udd(1,0,0)
    pprint_Gamma_udd(1,1,1)
    pprint_Gamma_udd(1,2,2)
    pprint_Gamma_udd(1,3,3)
    print
    pprint_Gamma_udd(2,2,1)
    pprint_Gamma_udd(2,1,2)
    pprint_Gamma_udd(2,3,3)
    print
    pprint_Gamma_udd(3,2,3)
    pprint_Gamma_udd(3,3,2)
    pprint_Gamma_udd(3,1,3)
    pprint_Gamma_udd(3,3,1)
    print"-"*40
    print"Ricci tensor:"
    pprint_Rmn_dd(0,0)
    e =  Rmn.dd(1,1)
    pprint_Rmn_dd(1,1)
    pprint_Rmn_dd(2,2)
    pprint_Rmn_dd(3,3)
    #print
    #print "scalar curvature:"
    #print curvature(Rmn)
    print "-"*40
    print "solve the Einstein's equations:"
    e = e.subs(nu(r), -lam(r))
    l =  dsolve(e, [lam(r)])
    pprint( Eq(lam(r), l) )
    metric = gdd.subs(lam(r), l).subs(nu(r),-l)#.combine()
    print "metric:"
    pprint( metric )

def derive_true_solution():
    import sympy as sym
    u = sym.symbols('u', cls=sym.Function)  # function u(t)
    t, a, p, I = sym.symbols('t a p I', real=True)

    def ode(u, t, a, p):
        """Define ODE: u' = (a + p*t)*u. Return residual."""
        return sym.diff(u, t) + (a + p*t)*u

    eq = ode(u(t), t, a, p)
    s = sym.dsolve(eq)
    # s is sym.Eq object u(t) == expression, we want u = expression,
    # so grab the right-hand side of the equality (Eq obj.)
    u = s.rhs
    print u
    # u contains C1, replace it with a symbol we can fit to
    # the initial condition
    C1 = sym.symbols('C1', real=True)
    u = u.subs('C1', C1)
    print u
    # Initial condition equation
    eq = u.subs(t, 0) - I
    s = sym.solve(eq, C1)  # solve eq wrt C1
    print s
    # s is a list s[0] = ...
    # Replace C1 in u by the solution
    u = u.subs(C1, s[0])
    print 'u:', u
    print sym.latex(u)  # latex formula for reports

    # Consistency check: u must fulfill ODE and initial condition
    print 'ODE is fulfilled:', sym.simplify(ode(u, t, a, p))
    print 'u(0)-I:', sym.simplify(u.subs(t, 0) - I)

    # Convert u expression to Python numerical function
    # (modules='numpy' allows numpy arrays as arguments,
    # we want this for t)
    u_func = sym.lambdify([t, I, a, p], u, modules='numpy')
    return u_func

def main():

    #print g
    print "-"*40
    print "Christoffel symbols:"
    print Gamma.udd(0,1,0)
    print Gamma.udd(0,0,1)
    print
    print Gamma.udd(1,0,0)
    print Gamma.udd(1,1,1)
    print Gamma.udd(1,2,2)
    print Gamma.udd(1,3,3)
    print
    print Gamma.udd(2,2,1)
    print Gamma.udd(2,1,2)
    print Gamma.udd(2,3,3)
    print
    print Gamma.udd(3,2,3)
    print Gamma.udd(3,3,2)
    print Gamma.udd(3,1,3)
    print Gamma.udd(3,3,1)
    print "-"*40
    print "Ricci tensor:"
    print Rmn.dd(0,0)
    e =  Rmn.dd(1,1)
    print e
    print Rmn.dd(2,2)
    print Rmn.dd(3,3)
    #print
    #print "scalar curvature:"
    #print curvature(Rmn)
    print "-"*40
    print "solve the Einstein's equations:"
    e = e.subs(nu(r), -lam(r))
    l =  dsolve(e, [lam(r)])
    print lam(r)," = ",l
    metric = gdd.subs(lam(r), l).subs(nu(r),-l)#.combine()
    print "metric:"
    print metric

    def __call__(self, equations, variables=None, method_options=None):
        logger.warn("The 'independent' state updater is deprecated and might be "
                    "removed in future versions of Brian.",
                    'deprecated_independent', once=True)
        method_options = extract_method_options(method_options, {})
        if equations.is_stochastic:
            raise UnsupportedEquationsException('Cannot solve stochastic '
                                                'equations with this state '
                                                'updater')
        if variables is None:
            variables = {}

        diff_eqs = equations.get_substituted_expressions(variables)

        t = Symbol('t', real=True, positive=True)
        dt = Symbol('dt', real=True, positive=True)
        t0 = Symbol('t0', real=True, positive=True)

        code = []
        for name, expression in diff_eqs:
            rhs = str_to_sympy(expression.code, variables)

            # We have to be careful and use the real=True assumption as well,
            # otherwise sympy doesn't consider the symbol a match to the content
            # of the equation
            var = Symbol(name, real=True)
            f = sp.Function(name)
            rhs = rhs.subs(var, f(t))
            derivative = sp.Derivative(f(t), t)
            diff_eq = sp.Eq(derivative, rhs)
            # TODO: simplify=True sometimes fails with 0.7.4, see:
            # https://github.com/sympy/sympy/issues/2666
            try:
                general_solution = sp.dsolve(diff_eq, f(t), simplify=True)
            except RuntimeError:
                general_solution = sp.dsolve(diff_eq, f(t), simplify=False)
            # Check whether this is an explicit solution
            if not getattr(general_solution, 'lhs', None) == f(t):
                raise UnsupportedEquationsException('Cannot explicitly solve: '
                                                    + str(diff_eq))
            # Solve for C1 (assuming "var" as the initial value and "t0" as time)
            if general_solution.has(Symbol('C1')):
                if general_solution.has(Symbol('C2')):
                    raise UnsupportedEquationsException('Too many constants in solution: %s' % str(general_solution))
                constant_solution = sp.solve(general_solution, Symbol('C1'))
                if len(constant_solution) != 1:
                    raise UnsupportedEquationsException(("Couldn't solve for the constant "
                                                         "C1 in : %s ") % str(general_solution))
                constant = constant_solution[0].subs(t, t0).subs(f(t0), var)
                solution = general_solution.rhs.subs('C1', constant)
            else:
                solution = general_solution.rhs.subs(t, t0).subs(f(t0), var)
            # Evaluate the expression for one timestep
            solution = solution.subs(t, t + dt).subs(t0, t)
            # only try symplifying it -- it sometimes raises an error
            try:
                solution = solution.simplify()
            except ValueError:
                pass

            code.append(name + ' = ' + sympy_to_str(solution))

        return '\n'.join(code)

def eq1():
    r = Symbol("r")
    e = relativity.Rmn.dd(0,0)
    e = e.subs(relativity.nu(r), -relativity.lam(r))
    print dsolve(e, [relativity.lam(r)])