def test_log_simplify():
    x = Symbol("x", positive=True)
    assert log(x**2).expand() == 2*log(x)
    assert expand_log(log(x**(2 + log(2)))) == (2 + log(2))*log(x)

    z = Symbol('z')
    assert log(sqrt(z)).expand() == log(z)/2
    assert expand_log(log(z**(log(2) - 1))) == (log(2) - 1)*log(z)
    assert log(z**(-1)).expand() != -log(z)
    assert log(z**(x/(x+1))).expand() == x*log(z)/(x + 1)

def test_issue_8866():
    assert simplify(log(x, 10, evaluate=False)) == simplify(log(x, 10))
    assert expand_log(log(x, 10, evaluate=False)) == expand_log(log(x, 10))

    y = Symbol('y', positive=True)
    l1 = log(exp(y), exp(10))
    b1 = log(exp(y), exp(5))
    l2 = log(exp(y), exp(10), evaluate=False)
    b2 = log(exp(y), exp(5), evaluate=False)
    assert simplify(log(l1, b1)) == simplify(log(l2, b2))
    assert expand_log(log(l1, b1)) == expand_log(log(l2, b2))

 def _eval_simplify(self, ratio, measure):
     from sympy.simplify.simplify import expand_log, simplify
     if (len(self.args) == 2):
         return simplify(self.func(*self.args), ratio=ratio, measure=measure)
     expr = self.func(simplify(self.args[0], ratio=ratio, measure=measure))
     expr = expand_log(expr, deep=True)
     return min([expr, self], key=measure)

 def _eval_simplify(self, ratio, measure, rational, inverse):
     from sympy.simplify.simplify import expand_log, simplify, inversecombine
     if (len(self.args) == 2):
         return simplify(self.func(*self.args), ratio=ratio, measure=measure,
                         rational=rational, inverse=inverse)
     expr = self.func(simplify(self.args[0], ratio=ratio, measure=measure,
                      rational=rational, inverse=inverse))
     if inverse:
         expr = inversecombine(expr)
     expr = expand_log(expr, deep=True)
     return min([expr, self], key=measure)

def canonise_log(equation):
    expanded_log = sympy.expand_log(equation, force=True)

    terms = expanded_log.as_ordered_terms()

    a, b = sympy.Wild('a'), sympy.Wild('b')

    total_interior = 1
    for term in terms:

        if sympy.ask(sympy.Q.complex(term)):
            term_interior *= -1
        else:
            term_interior = term.match(sympy.log(a) / b)[a]

            if term.could_extract_minus_sign():
                total_interior /= term_interior
            else:
                total_interior *= term_interior

    if isinstance(total_interior, sympy.Add):
        invert = False
    elif isinstance(total_interior, sympy.Mul):

        match = total_interior.together().match(x / b)  # for some reason, (x/3).match(a/b) gives {a: 1/3, b: 1/x} so we have to use a workaround
        if match is not None:
            invert = False
        else:
            match = total_interior.together().match(a / b)

            degree_numerator = 0 if isinstance(match[a], sympy.Rational) else match[a].as_poly().degree()
            degree_denominator = 0 if isinstance(match[b], sympy.Rational) else match[b].as_poly().degree()

            if degree_numerator < degree_denominator:
                invert = True
            else:
                invert = False
    elif isinstance(total_interior, sympy.Pow):
        index = total_interior.as_base_exp()[1]
        if index < 0:
            invert = True
        else:
            invert = False

    else:  # for debugging - wtf kind of c-c-c-class is it???
        print(total_interior, type(total_interior))

    if invert:
        return -sympy.log((1 / total_interior).together(), evaluate=False) / terms[0].as_coeff_Mul()[0].q
    else:
        return sympy.log(total_interior.together(), evaluate=False) / terms[0].as_coeff_Mul()[0].q

    def _eval_expand_log(self, deep=True, **hints):
        from sympy import unpolarify, expand_log
        from sympy.concrete import Sum, Product
        force = hints.get('force', False)
        if (len(self.args) == 2):
            return expand_log(self.func(*self.args), deep=deep, force=force)
        arg = self.args[0]
        if arg.is_Integer:
            # remove perfect powers
            p = perfect_power(int(arg))
            if p is not False:
                return p[1]*self.func(p[0])
        elif arg.is_Rational:
            return log(arg.p) - log(arg.q)
        elif arg.is_Mul:
            expr = []
            nonpos = []
            for x in arg.args:
                if force or x.is_positive or x.is_polar:
                    a = self.func(x)
                    if isinstance(a, log):
                        expr.append(self.func(x)._eval_expand_log(**hints))
                    else:
                        expr.append(a)
                elif x.is_negative:
                    a = self.func(-x)
                    expr.append(a)
                    nonpos.append(S.NegativeOne)
                else:
                    nonpos.append(x)
            return Add(*expr) + log(Mul(*nonpos))
        elif arg.is_Pow or isinstance(arg, exp):
            if force or (arg.exp.is_real and (arg.base.is_positive or ((arg.exp+1)
                .is_positive and (arg.exp-1).is_nonpositive))) or arg.base.is_polar:
                b = arg.base
                e = arg.exp
                a = self.func(b)
                if isinstance(a, log):
                    return unpolarify(e) * a._eval_expand_log(**hints)
                else:
                    return unpolarify(e) * a
        elif isinstance(arg, Product):
            if arg.function.is_positive:
                return Sum(log(arg.function), *arg.limits)

        return self.func(arg)

def test_risch_integrate():
    assert risch_integrate(t0*exp(x), x) == t0*exp(x)
    assert risch_integrate(sin(x), x, rewrite_complex=True) == -exp(I*x)/2 - exp(-I*x)/2

    # From my GSoC writeup
    assert risch_integrate((1 + 2*x**2 + x**4 + 2*x**3*exp(2*x**2))/
    (x**4*exp(x**2) + 2*x**2*exp(x**2) + exp(x**2)), x) == \
        NonElementaryIntegral(exp(-x**2), x) + exp(x**2)/(1 + x**2)


    assert risch_integrate(0, x) == 0

    # also tests prde_cancel()
    e1 = log(x/exp(x) + 1)
    ans1 = risch_integrate(e1, x)
    assert ans1 == (x*log(x*exp(-x) + 1) + NonElementaryIntegral((x**2 - x)/(x + exp(x)), x))
    assert cancel(diff(ans1, x) - e1) == 0

    # also tests issue #10798
    e2 = (log(-1/y)/2 - log(1/y)/2)/y - (log(1 - 1/y)/2 - log(1 + 1/y)/2)/y
    ans2 = risch_integrate(e2, y)
    assert ans2 == log(1/y)*log(1 - 1/y)/2 - log(1/y)*log(1 + 1/y)/2 + \
            NonElementaryIntegral((I*pi*y**2 - 2*y*log(1/y) - I*pi)/(2*y**3 - 2*y), y)
    assert expand_log(cancel(diff(ans2, y) - e2), force=True) == 0

    # These are tested here in addition to in test_DifferentialExtension above
    # (symlogs) to test that backsubs works correctly.  The integrals should be
    # written in terms of the original logarithms in the integrands.

    # XXX: Unfortunately, making backsubs work on this one is a little
    # trickier, because x**x is converted to exp(x*log(x)), and so log(x**x)
    # is converted to x*log(x). (x**2*log(x)).subs(x*log(x), log(x**x)) is
    # smart enough, the issue is that these splits happen at different places
    # in the algorithm.  Maybe a heuristic is in order
    assert risch_integrate(log(x**x), x) == x**2*log(x)/2 - x**2/4

    assert risch_integrate(log(x**y), x) == x*log(x**y) - x*y
    assert risch_integrate(log(sqrt(x)), x) == x*log(sqrt(x)) - x/2

def test_log_simplify():
    x = Symbol("x", positive=True)
    assert log(x**2).expand() == 2*log(x)
    assert expand_log(log(x**(2+log(2)))) == (2+log(2))*log(x)

                new_exp_terms.append(exp_term)
            else:
                done = True
                new_exp_terms.append(attempt)
    if not done:
        new_exp_terms.append(numsum)
    return e.func(*chain(new_exp_terms, non_num_other))


expm1_opt = ReplaceOptim(lambda e: e.is_Add, _expm1_value)


log1p_opt = ReplaceOptim(
    lambda e: isinstance(e, log),
    lambda l: expand_log(l.replace(
        log, lambda arg: log(arg.factor())
    )).replace(log(_u+1), log1p(_u))
)

def create_expand_pow_optimization(limit):
    """ Creates an instance of :class:`ReplaceOptim` for expanding ``Pow``.

    The requirements for expansions are that the base needs to be a symbol
    and the exponent needs to be an integer (and be less than or equal to
    ``limit``).

    Parameters
    ==========

    limit : int
         The highest power which is expanded into multiplication.

# <codecell>

import sympy
from sympy.abc import x, z
p=sympy.symbols('p',positive=True)
sympy.init_printing()

# <codecell>

L=p**x*(1-p)**(1-x)
J=np.prod([L.subs(x,i) for i in xs]) # objective function to maximize
J

# <codecell>

logJ=sympy.expand_log(sympy.log(J))
sol=sympy.solve(sympy.diff(logJ,p),p)[0]

x=linspace(0,1,100)
plot(x,map(sympy.lambdify(p,logJ,'numpy'),x),sol,logJ.subs(p,sol),'o',
                                          p_true,logJ.subs(p,p_true),'s',)
xlabel('$p$',fontsize=18)
ylabel('Likelihood',fontsize=18)
title('Estimate not equal to true value',fontsize=18)

# <codecell>

L=p**x*(1-p)**(1-x)
J=np.prod([L.subs(x,i) for i in xs]) # objective function