def test_together2():
    x, y, z = symbols("xyz")
    assert together(1/(x*y) + 1/y**2) == 1/x*y**(-2)*(x + y)
    assert together(1/(1 + 1/x)) == x/(1 + x)
    x = symbols("x", nonnegative=True)
    y = symbols("y", real=True)
    assert together(1/x**y + 1/x**(y-1)) == x**(-y)*(1 + x)

def showAnswer(p):
    n = Symbol('n')
    sp = together(sumFibPower(n, p, viaPowers=True).expand().simplify())
    sn = together(sumFibPower(n, p, viaPowers=False).expand().simplify())

    print ("Sum F(i)^%d = "%(p))
    print ("=",sp,"=")
    print ("=",sn)

def singular_term(F,X,Y,L,I,version):
    """
    Computes a single set of singular terms of the Puiseux expansions.

    For `I=1`, the function computes the first term of each finite 
    `\mathbb{K}` expansion. For `I=2`, it computes the other terms of the
    expansion.
    """
    T = []

    # if the curve is singular then compute the singular tuples
    # otherwise, use the standard newton polygon method
    if is_singular(F,X,Y):
        for (q,m,l,Phi) in desingularize(F,X,Y):
            for eta in Phi.all_roots(radicals=True):
                tau = (q,1,m,1,sympy.together(eta))
                T.append((tau,0,1))
    else:
        # each side of the newton polygon corresponds to a K-term.
        for (q,m,l,Phi) in polygon(F,X,Y,I):
            # the rational method
            if version == 'rational':
                u,v = _bezout(q,m)      

            # each newton polygon side has a characteristic
            # polynomial. For each square-free factor, each root
            # corresponds to a K-term
            for (Psi,r) in _square_free(Phi,_Z):
                Psi = sympy.Poly(Psi,_Z)
                
                # compute the roots of Psi. Use the RootOf construct if 
                # possible. In the case when Psi is over EX (i.e. when
                # RootOf doesn't work) then compute symbolic roots.
                try:
                    roots = Psi.all_roots(radicals=True)
                except NotImplementedError:
                    roots = sympy.roots(Psi,_Z).keys()

                for xi in roots:
                    # the classical version returns the "raw" roots
                    if version == 'classical':
                        P = sympy.Poly(_U**q-xi,_U)
                        beta = RootOf(P,0,radicals=True)
                        tau = (q,1,m,beta,1)
                        T.append((tau,l,r))
                    # the rational version rescales parameters so as
                    # to include ony rational terms in the Puiseux
                    # expansions.
                    if version == 'rational':
                        mu = xi**(-v)
                        beta = sympy.together(xi**u)
                        tau = (q,mu,m,beta,1)
                        T.append((tau,l,r))
    return T

 def __mul__(self, other):
     if isinstance(other, RayTransferMatrix):
         return RayTransferMatrix(Matrix.__mul__(self, other))
     elif isinstance(other, GeometricRay):
         return GeometricRay(Matrix.__mul__(self, other))
     elif isinstance(other, BeamParameter):
         temp = self*Matrix(((other.q,), (1,)))
         q = (temp[0]/temp[1]).expand(complex=True)
         return BeamParameter(other.wavelen,
                              together(re(q)),
                              z_r=together(im(q)))
     else:
         return Matrix.__mul__(self, other)

def findSymbolicExpressionL1(df,target,result_object,scaler,task="regression",logged_target=False,find_acceptable=True,features_type='best_features',max_length=-1):
    columns=df.columns.values
    X=result_object[features_type]
    individuals=result_object['best_individuals_equations']

    
    final=[]
    eqs=convertIndividualsToEqs(individuals,columns)
    eqs2=[]
    if max_length>0:
        for eq in eqs:
            if len(eq)<max_length:
                eqs2.append(eq)           
        eqs=eqs2
    
    x=np.column_stack(X)        
    if task=='regression':
        models=findSymbolicExpressionL1_regression_helper(x,target)
    elif task=='classification':
        models=findSymbolicExpressionL1_classification_helper(x,target)  
        
    
    
    
    for model,score,coefs in models:
        if(len(model.coef_.shape)==1):                              
            eq=""
            #category represents a class. In logistic regression with scikit learn                
            for j in range(0,len(eqs)):
                    eq=eq+"+("+str(coefs[j])+"*("+str(sympify(eqs[j]))+"))"
                                  
            final.append((together(sympify(eq)),np.around(score,3)))
        else:
            dummy=[]
            for category_coefs in model.coef_:               
                eq=""
                #category represents a class. In logistic regression with scikit learn
                
                for j in range(0,len(eqs)):    
                    eq=eq+"+("+str(category_coefs[j])+"*("+str(sympify(eqs[j]))+"))"
                    
                dummy.append(together(sympify(eq)))
            final.append((dummy,np.around(score,3)))
    
    
    if find_acceptable:
        final=findAcceptableL1(final)
    
    return final     

def antal_h_coefficient(index, game_matrix):
    """
    Returns the H_index coefficient, according to Antal et al. (2009), as given by equation 2.
    H_k = \frac{1}{n^2} \sum_{i=1}^{n} \sum_{j=1}^{n} (a_{kj}-a_{jj})
    Parameters
    ----------
    index: int
    game_matrix: sympy.Matrix
    
    Returns
    -------
    out: sympy.Expr
    
    
    Examples:
    --------
    >>>a = Symbol('a')
    >>>antal_h_coefficient(0, Matrix([[a,2,3],[4,5,6],[7,8,9]]))
    Out[1]: (2*a - 29)/9
    
    >>>antal_h_coefficient(0, Matrix([[1,2,3],[4,5,6],[7,8,9]]))
    Out[1]: -3
    """
    size = game_matrix.shape[0]
    suma = 0
    for i in range(0, size):
        for j in range(0, size):
            suma = suma + (game_matrix[index, i] - game_matrix[i, j])
    return sympy.together(sympy.simplify(suma / (size ** 2)), size)

def _simplify(expr):
    u"""Simplifie une expression.

    Alias de simplify (sympy 0.6.4).
    Mais simplify n'est pas garanti d'être stable dans le temps.
    (Cf. simplify.__doc__)."""
    return together(expand(Poly.cancel(powsimp(expr))))

def antal_l_coefficient(index, game_matrix):
    """
    Returns the L_index coefficient, according to Antal et al. (2009), as given by equation 1.
    L_k = \frac{1}{n} \sum_{i=1}^{n} (a_{kk}+a_{ki}-a_{ik}-a_{ii})
    Parameters
    ----------
    index: int
    game_matrix: sympy.Matrix
    
    Returns
    -------
    out: sympy.Expr
    
    
    Examples:
    --------
    >>>a = Symbol('a')
    >>>antal_l_coefficient(0, Matrix([[a,2,3],[4,5,6],[7,8,9]]))
    Out[1]: 2*(a - 10)/3
    
    
    >>>antal_l_coefficient(0, Matrix([[1,2,3],[4,5,6],[7,8,9]]))
    Out[1]: -6
    """
    size = game_matrix.shape[0]
    suma = 0
    for i in range(0, size):
        suma = suma + (game_matrix[index, index] + game_matrix[index, i] - game_matrix[i, index] - game_matrix[i, i])
    return sympy.together(sympy.simplify(suma / size), size)

def calculate_gains(sol, xin, optimize=True):
    """Calculate low-frequency gain and roots of a transfer function.

    **Parameters:**

    sol : dict
        the circuit solution
    xin : Sympy symbol
        the input variable
    optimize : boolean, optional
        If ``optimize`` is set to ``False``, no algebraic simplification
        will be attempted on the results. The default (``optimize=True``)
        results in ``sympy.together`` being called on each expression.

    **Returns:**

    gs : dict
        A dictionary with as keys the strings <key>/<xin> and as values
        dictionaries with keys ``'gain'``, ``'gain0'``, ``'poles'``,
        ``'zeros'``.

    """
    gains = {}
    for key, value in sol.items():
        tf = {}
        gain = sympy.together(value.diff(xin)) if optimize else value.diff(xin)
        (ps, zs) = get_roots(gain)
        tf.update({'gain': gain})
        tf.update({'gain0': gain.subs(s, 0)})
        tf.update({'poles': ps})
        tf.update({'zeros': zs})
        gains.update({"%s/%s" % (str(key), str(xin)): tf})
    return gains

 def apply(self, expr, evaluation):
     'Together[expr_]'
     
     expr_sympy = expr.to_sympy()
     result = sympy.together(expr_sympy)
     result = from_sympy(result)
     result = cancel(result)
     return result

def regular(S,X,Y,nterms,degree_bound):
    """
    INPUT:

    -- ``S``: a finite set of pairs `\{(\pi_k,F_k)\}_{1 \leq k \leq B}` 
         where `\pi_k` is a finite `\mathbb{K}`-expansion, `F_k
         \in \bar{\mathbb{K}}[X,Y]` with `F_k(0,0) = 0`, `\partial_Y
         F_k(0,0) \neq 0`, and `F_k(X,0) \neq 0`.

    -- ``H``: a positive integer

    OUTPUT:
    
    -- ``(list)``: a set `\{ \pi_k' \}_{1 \leq k \leq B}` of finite
         `\mathbb{K}` expansions such that `\pi_k'` begins with
         `\pi_k` and contains at least `H` `\mathbb{K}`-terms.

    """
    R = []
    for (pi,F) in S:        
        # grow each expansion to the number of desired terms
        Npi = len(pi)

        # if a degree bound is specified, get the degree of the
        # singular part of the Puiseux series. We also want to compute
        # enough terms to distinguish the puiseux series.
        q,mu,m,beta,eta = pi[-1]
        e = reduce( lambda q1,q2: q1*q1, (tau[0] for tau in pi) )
        ydeg = sum( tau[0]*tau[2] for tau in pi )
        need_more_terms = True

        while ((Npi < nterms) and (ydeg < e*degree_bound)) or need_more_terms:
            # if the set of all (0,j), j!=0 is empty, then we've 
            # encountered a finite puiseux expansion
            a = dict(F.terms())
            ms = [j for (j,i) in a.keys() if i==0 and j!=0]
            if ms == []: break
            else:        m = min(ms)

            # if a degree bound is specified, break pi-series
            # construction once we break the bound.
            ydeg += m
            Npi += 1
                
            beta = sympy.together(-a[(m,0)]/a[(0,1)])
            tau = (1,1,m,beta,1)
            pi.append(tau)
            F = _new_polynomial(F,X,Y,tau,m)

            # if the degree in the y-series isn't divisible by the
            # ramification index then we have enough terms to
            # distinguish the puiseux series.
            if sympy.gcd(ydeg,e) == 1:
                need_more_terms = False

        R.append(tuple(pi))

    return tuple(R)

def test_apart_extension():
    f = 2/(x**2 + 1)
    g = I/(x + I) - I/(x - I)

    assert apart(f, extension=I) == g
    assert apart(f, gaussian=True) == g

    f = x/((x - 2)*(x + I))

    assert factor(together(apart(f))) == f

def exercise_6_43():

    z = sp.Symbol('z')

    b_of_z = 1/(1 - z)

    sp.pprint('\n' + 'b_of_z:')
    sp.pprint(b_of_z)

    e_b_prime_of_z = z**2/(1 - z)**2

    sp.pprint('\n' +'e_b_prime_of_z:')
    sp.pprint(e_b_prime_of_z)

    e_b_prime_of_n = sp.series(e_b_prime_of_z, z, 0, 10)

    sp.pprint('\n' +'e_b_prime_of_n:')
    sp.pprint(e_b_prime_of_n)

    e_b_of_z = sp.integrate(e_b_prime_of_z, z)
    constant = e_b_of_z.subs(z, 0)
    e_b_of_z = e_b_of_z - constant

    sp.pprint('\n' +'e_b_of_z:')
    sp.pprint(e_b_of_z)

    e_b_of_n = sp.series(e_b_of_z, z, 0, 10)

    sp.pprint('\n' +'e_b_of_n:')
    sp.pprint(e_b_of_n)

    e_b_prime_of_z = sp.diff(e_b_of_z, z)

    sp.pprint('\n' +'e_b_prime_of_z:')
    sp.pprint(e_b_prime_of_z)

    f_of_z = sp.together((1 - z)**2*e_b_prime_of_z)

    sp.pprint('\n' +'f_of_z:')
    sp.pprint(f_of_z)

    F_of_z = sp.integrate(f_of_z, z)

    sp.pprint('\n' +'F_of_z:')
    sp.pprint(F_of_z)

    c_b_of_z = (e_b_of_z.subs(z, 0) + F_of_z)/(1 - z)**2

    sp.pprint('\n' +'c_b_of_z:')
    sp.pprint(c_b_of_z)

    c_b_of_n = sp.series(c_b_of_z, z, 0, 10)

    sp.pprint('\n' +'c_b_of_n:')
    sp.pprint(c_b_of_n)

def findSymbolicExpression(df,target,result_object,scaler,task="regression",logged_target=False,acceptable_only=True):
    columns=df.columns.values
    X=result_object['best_features']
    individuals=result_object['best_individuals_object']

    
    final=[]
    eqs=convertIndividualsToEqs(individuals,columns)
    for i in range(1,len(X)):
        x=np.column_stack(X[0:i])        
        if task=='regression':
            model=linear_model.LinearRegression()
            scoring='r2'
        elif task=='classification':
            model=linear_model.LogisticRegression()  
            scoring='f1_weighted'
        
        scores=cross_validation.cross_val_score(estimator=model, X=x, y=target, cv=10,scoring=scoring)
        model.fit(x,target)  
        
        m=np.mean(scores)
        sd=np.std(scores)
        if(acceptable_only and m>sd) or (not acceptable_only):
            if(len(model.coef_.shape)==1):                              
                eq=""
                #category represents a class. In logistic regression with scikit learn                
                for j in range(0,i):
                        eq=eq+"+("+str(model.coef_[j])+"*("+str(sympify(eqs[j]))+"))"
                                    
                final.append((together(sympify(eq)),m,sd))
            else:
                dummy=[]
                for category_coefs in model.coef_:               
                    eq=""
                    #category represents a class. In logistic regression with scikit learn
                    
                    for j in range(0,i):    
                        eq=eq+"+("+str(category_coefs[j])+"*("+str(sympify(eqs[j]))+"))"
                        
                    dummy.append(together(sympify(eq)))
                final.append((dummy,m,sd))
    return final   

def sympy_factor(expr_sympy):
    try:
        result = sympy.together(expr_sympy)
        numer, denom = result.as_numer_denom()
        if denom == 1:
            result = sympy.factor(expr_sympy)
        else:
            result = sympy.factor(numer) / sympy.factor(denom)
    except sympy.PolynomialError:
        return expr_sympy
    return result

 def apply(self, expr, evaluation):
     'Simplify[expr_]'
     
     expr_sympy = expr.to_sympy()
     result = expr_sympy
     result = sympy.simplify(result)
     result = sympy.trigsimp(result)
     result = sympy.together(result)
     result = sympy.cancel(result)
     result = from_sympy(result)
     return result

def calculate_gains(sol, xin, optimize=True):
    gains = {}
    for key, value in sol.iteritems():
        tf = {}
        gain = sympy.together(value.diff(xin)) if optimize else value.diff(xin)
        (ps, zs) = get_roots(gain)
        tf.update({'gain': gain})
        tf.update({'gain0': gain.subs(s, 0)})
        tf.update({'poles': ps})
        tf.update({'zeros': zs})
        gains.update({"%s/%s" % (str(key), str(xin)): tf})
    return gains

    def _get_Ylm(self, l, m):
        """
        Compute an expression for spherical harmonic of order (l,m)
        in terms of Cartesian unit vectors, :math:`\hat{z}`
        and :math:`\hat{x} + i \hat{y}`

        Parameters
        ----------
        l : int
            the degree of the harmonic
        m : int
            the order of the harmonic; |m| < l

        Returns
        -------
        expr :
            a sympy expression that corresponds to the
            requested Ylm

        References
        ----------
        https://en.wikipedia.org/wiki/Spherical_harmonics
        """
        import sympy as sp

        # the relevant cartesian and spherical symbols
        x, y, z, r = sp.symbols('x y z r', real=True, positive=True)
        xhat, yhat, zhat = sp.symbols('xhat yhat zhat', real=True, positive=True)
        xpyhat = sp.Symbol('xpyhat', complex=True)
        phi, theta = sp.symbols('phi theta')
        defs = [(sp.sin(phi), y/sp.sqrt(x**2+y**2)),
                (sp.cos(phi), x/sp.sqrt(x**2+y**2)),
                (sp.cos(theta), z/sp.sqrt(x**2 + y**2 + z**2))
                ]

        # the cos(theta) dependence encoded by the associated Legendre poly
        expr = sp.assoc_legendre(l, m, sp.cos(theta))

        # the exp(i*m*phi) dependence
        expr *= sp.expand_trig(sp.cos(m*phi)) + sp.I*sp.expand_trig(sp.sin(m*phi))

        # simplifying optimizations
        expr = sp.together(expr.subs(defs)).subs(x**2 + y**2 + z**2, r**2)
        expr = expr.expand().subs([(x/r, xhat), (y/r, yhat), (z/r, zhat)])
        expr = expr.factor().factor(extension=[sp.I]).subs(xhat+sp.I*yhat, xpyhat)
        expr = expr.subs(xhat**2 + yhat**2, 1-zhat**2).factor()

        # and finally add the normalization
        amp = sp.sqrt((2*l+1) / (4*numpy.pi) * sp.factorial(l-m) / sp.factorial(l+m))
        expr *= amp

        return expr

def test_action_verbs():
    assert nsimplify((1/(exp(3*pi*x/5)+1))) == (1/(exp(3*pi*x/5)+1)).nsimplify()
    assert ratsimp(1/x + 1/y) == (1/x + 1/y).ratsimp()
    assert trigsimp(log(x), deep=True) == (log(x)).trigsimp(deep = True)
    assert radsimp(1/(2+sqrt(2))) == (1/(2+sqrt(2))).radsimp()
    assert powsimp(x**y*x**z*y**z, combine='all') == (x**y*x**z*y**z).powsimp(combine='all')
    assert simplify(x**y*x**z*y**z) == (x**y*x**z*y**z).simplify()
    assert together(1/x + 1/y) == (1/x + 1/y).together()
    assert separate((x*(y*z)**3)**2) == ((x*(y*z)**3)**2).separate()
    assert collect(a*x**2 + b*x**2 + a*x - b*x + c, x) == (a*x**2 + b*x**2 + a*x - b*x + c).collect(x)
    assert apart(y/(y+2)/(y+1), y) == (y/(y+2)/(y+1)).apart(y)
    assert combsimp(y/(x+2)/(x+1)) == (y/(x+2)/(x+1)).combsimp()
    assert factor(x**2+5*x+6) == (x**2+5*x+6).factor()
    assert refine(sqrt(x**2)) == sqrt(x**2).refine()
    assert cancel((x**2+5*x+6)/(x+2)) == ((x**2+5*x+6)/(x+2)).cancel()

 def apply_list(self, expr, vars, evaluation):
     'FactorTermsList[expr_, vars_List]'
     if expr == Integer(0):
         return Expression('List', Integer(1), Integer(0))
     elif isinstance(expr, Number):
         return Expression('List', expr, Integer(1))
     
     for x in vars.leaves:
         if not(isinstance(x, Atom)):
             return evaluation.message('CoefficientList', 'ivar', x)
     
     sympy_expr = expr.to_sympy()
     if sympy_expr is None:
         return Expression('List', Integer(1), expr)
     sympy_expr = sympy.together(sympy_expr)
     
     sympy_vars = [x.to_sympy() for x in vars.leaves if isinstance(x, Symbol) and sympy_expr.is_polynomial(x.to_sympy())]
     
     result = []
     numer, denom = sympy_expr.as_numer_denom()
     try:
         from sympy import factor, factor_list, Poly
         if denom == 1:
             # Get numerical part
             num_coeff, num_polys = factor_list(Poly(numer))
             result.append(num_coeff)
             
             # Get factors are independent of sub list of variables
             if (sympy_vars and isinstance(expr, Expression) 
                 and any(x.free_symbols.issubset(sympy_expr.free_symbols) for x in sympy_vars)):
                 for i in reversed(range(len(sympy_vars))):
                     numer = factor(numer) / factor(num_coeff)
                     num_coeff, num_polys = factor_list(Poly(numer), *[x for x in sympy_vars[:(i+1)]])
                     result.append(sympy.expand(num_coeff))
             
             # Last factor
             numer = factor(numer) / factor(num_coeff)
             result.append(sympy.expand(numer))
         else:
             num_coeff, num_polys = factor_list(Poly(numer))
             den_coeff, den_polys = factor_list(Poly(denom))
             result = [num_coeff / den_coeff, sympy.expand(factor(numer)/num_coeff / (factor(denom)/den_coeff))]
     except sympy.PolynomialError: # MMA does not raise error for non poly
         result.append(sympy.expand(numer))
         # evaluation.message(self.get_name(), 'poly', expr)
     
     return Expression('List', *[from_sympy(i) for i in result])