python类Eq()的实例源码

def error_bound(n):
    e = Symbol('e')
    y = Symbol('y')
    growth_function_bound = generate_growth_function_bound(50)
    a = original_vc_bound(n, 0.05, growth_function_bound)
    b = nsolve(Eq(rademacher_penalty_bound(n, 0.05, growth_function_bound), y), 1)
    c = nsolve(Eq(parrondo_van_den_broek_right(e, n, 0.05, growth_function_bound), e), 1)
    d = nsolve(Eq(devroye(e, n, 0.05, growth_function_bound), e), 1)
    return a, b, c, d


# def test_three():
#     e = Symbol('e')
#     y = Symbol('y')
#     n = Symbol('n')
#     growth_function_bound = generate_growth_function_bound(50)
#     a = original_vc_bound(5, 0.05, growth_function_bound)
#     b = nsolve(Eq(rademacher_penalty_bound(5, 0.05, growth_function_bound), y), 5)
#     c = nsolve(Eq(parrondo_van_den_broek_right(e, 5, 0.05, growth_function_bound), e), 1)
#     d = nsolve(Eq(devroye(e, 5, 0.05, growth_function_bound), e), 1)
#     return a, b, c, d

def test_simplify_other():
    assert simplify(sin(x)**2 + cos(x)**2) == 1
    assert simplify(gamma(x + 1)/gamma(x)) == x
    assert simplify(sin(x)**2 + cos(x)**2 + factorial(x)/gamma(x)) == 1 + x
    assert simplify(
        Eq(sin(x)**2 + cos(x)**2, factorial(x)/gamma(x))) == Eq(1, x)
    nc = symbols('nc', commutative=False)
    assert simplify(x + x*nc) == x*(1 + nc)
    # issue 6123
    # f = exp(-I*(k*sqrt(t) + x/(2*sqrt(t)))**2)
    # ans = integrate(f, (k, -oo, oo), conds='none')
    ans = I*(-pi*x*exp(-3*I*pi/4 + I*x**2/(4*t))*erf(x*exp(-3*I*pi/4)/
        (2*sqrt(t)))/(2*sqrt(t)) + pi*x*exp(-3*I*pi/4 + I*x**2/(4*t))/
        (2*sqrt(t)))*exp(-I*x**2/(4*t))/(sqrt(pi)*x) - I*sqrt(pi) * \
        (-erf(x*exp(I*pi/4)/(2*sqrt(t))) + 1)*exp(I*pi/4)/(2*sqrt(t))
    assert simplify(ans) == -(-1)**(S(3)/4)*sqrt(pi)/sqrt(t)
    # issue 6370
    assert simplify(2**(2 + x)/4) == 2**x

def __init__(self, expr, var_start_end_x, var_start_end_y,
            has_equality, use_interval_math, depth, nb_of_points):
        super(ImplicitSeries, self).__init__()
        self.expr = sympify(expr)
        self.var_x = sympify(var_start_end_x[0])
        self.start_x = float(var_start_end_x[1])
        self.end_x = float(var_start_end_x[2])
        self.var_y = sympify(var_start_end_y[0])
        self.start_y = float(var_start_end_y[1])
        self.end_y = float(var_start_end_y[2])
        self.get_points = self.get_raster
        self.has_equality = has_equality  # If the expression has equality, i.e.
                                         #Eq, Greaterthan, LessThan.
        self.nb_of_points = nb_of_points
        self.use_interval_math = use_interval_math
        self.depth = 4 + depth

def power_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if symbol not in exp.free_symbols and isinstance(base, sympy.Symbol):
        if sympy.simplify(exp + 1) == 0:
            return ReciprocalRule(base, integrand, symbol)
        return PowerRule(base, exp, integrand, symbol)
    elif symbol not in base.free_symbols and isinstance(exp, sympy.Symbol):
        rule = ExpRule(base, exp, integrand, symbol)

        if sympy.ask(~sympy.Q.zero(sympy.log(base))):
            return rule
        elif sympy.ask(sympy.Q.zero(sympy.log(base))):
            return ConstantRule(1, 1, symbol)

        return PiecewiseRule([
            (ConstantRule(1, 1, symbol), sympy.Eq(sympy.log(base), 0)),
            (rule, True)
        ], integrand, symbol)

def test_threaded():
    @threaded
    def function(expr, *args):
        return 2*expr + sum(args)

    assert function(Matrix([[x, y], [1, x]]), 1, 2) == \
        Matrix([[2*x + 3, 2*y + 3], [5, 2*x + 3]])

    assert function(Eq(x, y), 1, 2) == Eq(2*x + 3, 2*y + 3)

    assert function([x, y], 1, 2) == [2*x + 3, 2*y + 3]
    assert function((x, y), 1, 2) == (2*x + 3, 2*y + 3)

    assert function(set([x, y]), 1, 2) == set([2*x + 3, 2*y + 3])

    @threaded
    def function(expr, n):
        return expr**n

    assert function(x + y, 2) == x**2 + y**2
    assert function(x, 2) == x**2

def _eval_rewrite_as_Probability(self, arg, condition=None):
        rvs = arg.atoms(RandomSymbol)
        if len(rvs) > 1:
            raise NotImplementedError()
        if len(rvs) == 0:
            return arg

        rv = rvs.pop()
        if rv.pspace is None:
            raise ValueError("Probability space not known")

        symbol = rv.symbol
        if symbol.name[0].isupper():
            symbol = Symbol(symbol.name.lower())
        else :
            symbol = Symbol(symbol.name + "_1")

        if rv.pspace.is_Continuous:
            return Integral(arg.replace(rv, symbol)*Probability(Eq(rv, symbol), condition), (symbol, rv.pspace.domain.set.inf, rv.pspace.domain.set.sup))
        else:
            if rv.pspace.is_Finite:
                raise NotImplemented
            else:
                return Sum(arg.replace(rv, symbol)*Probability(Eq(rv, symbol), condition), (symbol, rv.pspace.domain.set.inf, rv.pspace.set.sup))

def test_euler_high_order():
    # an example from hep-th/0309038
    m = Symbol('m')
    k = Symbol('k')
    x = Function('x')
    y = Function('y')
    t = Symbol('t')
    L = (m*D(x(t), t)**2/2 + m*D(y(t), t)**2/2 -
         k*D(x(t), t)*D(y(t), t, t) + k*D(y(t), t)*D(x(t), t, t))
    assert euler(L, [x(t), y(t)]) == [Eq(2*k*D(y(t), t, t, t) -
                                         m*D(x(t), t, t)),
                                      Eq(-2*k*D(x(t), t, t, t) -
                                         m*D(y(t), t, t))]

    w = Symbol('w')
    L = D(x(t, w), t, w)**2/2
    assert euler(L) == [Eq(D(x(t, w), t, t, w, w))]

def __init__(self, expr, var_start_end_x, var_start_end_y,
            has_equality, use_interval_math, depth, nb_of_points,
            line_color):
        super(ImplicitSeries, self).__init__()
        self.expr = sympify(expr)
        self.var_x = sympify(var_start_end_x[0])
        self.start_x = float(var_start_end_x[1])
        self.end_x = float(var_start_end_x[2])
        self.var_y = sympify(var_start_end_y[0])
        self.start_y = float(var_start_end_y[1])
        self.end_y = float(var_start_end_y[2])
        self.get_points = self.get_raster
        self.has_equality = has_equality  # If the expression has equality, i.e.
                                         #Eq, Greaterthan, LessThan.
        self.nb_of_points = nb_of_points
        self.use_interval_math = use_interval_math
        self.depth = 4 + depth
        self.line_color = line_color

def power_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if symbol not in exp.free_symbols and isinstance(base, sympy.Symbol):
        if sympy.simplify(exp + 1) == 0:
            return ReciprocalRule(base, integrand, symbol)
        return PowerRule(base, exp, integrand, symbol)
    elif symbol not in base.free_symbols and isinstance(exp, sympy.Symbol):
        rule = ExpRule(base, exp, integrand, symbol)

        if fuzzy_not(sympy.log(base).is_zero):
            return rule
        elif sympy.log(base).is_zero:
            return ConstantRule(1, 1, symbol)

        return PiecewiseRule([
            (ConstantRule(1, 1, symbol), sympy.Eq(sympy.log(base), 0)),
            (rule, True)
        ], integrand, symbol)

def test_threaded():
    @threaded
    def function(expr, *args):
        return 2*expr + sum(args)

    assert function(Matrix([[x, y], [1, x]]), 1, 2) == \
        Matrix([[2*x + 3, 2*y + 3], [5, 2*x + 3]])

    assert function(Eq(x, y), 1, 2) == Eq(2*x + 3, 2*y + 3)

    assert function([x, y], 1, 2) == [2*x + 3, 2*y + 3]
    assert function((x, y), 1, 2) == (2*x + 3, 2*y + 3)

    assert function({x, y}, 1, 2) == {2*x + 3, 2*y + 3}

    @threaded
    def function(expr, n):
        return expr**n

    assert function(x + y, 2) == x**2 + y**2
    assert function(x, 2) == x**2

def _stamp_port(n, port, context):
    """
    Stamp a port in the MNA matrix by adding KVL and KCL equations to the
    kvl_equations and kcl_equations lists in the context dictionary.
    """
    an = context["a"][n]
    bn = context["b"][n]
    node_voltages = context["node_voltages"]
    kvl_equations = context["kvl_equations"]
    kcl_equations = context["kcl_equations"]

    node1, node2 = port["nodes"]
    v1 = node_voltages[node1]
    v2 = node_voltages[node2]
    port_voltage = (an + bn) / 2
    kvl_equation = sympy.Eq(port_voltage, v2 - v1)
    kvl_equations.append(kvl_equation)

    port_current = (an - bn) / (2 * context["port_resistances"][n])
    context["kcl_equations"][node1] += -port_current
    context["kcl_equations"][node2] += port_current

def _stamp_vcvs(n, vcvs, context):
    """
    Stamp a VCVS in the MNA matrix by adding KVL and KCL equations to the
    kvl_equations and kcl_equations lists in the context dictionary.
    """
    node_voltages = context["node_voltages"]
    kvl_equations = context["kvl_equations"]
    kcl_equations = context["kcl_equations"]
    vcvs_currents = context["vcvs_currents"]

    node1, node2, node3, node4 = vcvs["nodes"]
    v1 = node_voltages[node1]
    v2 = node_voltages[node2]
    v3 = node_voltages[node3]
    v4 = node_voltages[node4]
    gain = vcvs.get("gain", sympy.Symbol("A"))

    kvl_equation = sympy.Eq(gain * (v2 - v1), v4 - v3)
    kvl_equations.append(kvl_equation)

    current = vcvs_currents[n]
    kcl_equations[node3] += -current
    kcl_equations[node4] += current

def q_affine(expr, vars):
    """
    Return True if ``expr`` is (separately) affine in the variables ``vars``,
    False otherwise.

    Readapted from: https://stackoverflow.com/questions/36283548\
        /check-if-an-equation-is-linear-for-a-specific-set-of-variables/
    """
    # A function is (separately) affine in a given set of variables if all
    # non-mixed second order derivatives are identically zero.
    for x in as_tuple(vars):
        if x not in expr.atoms():
            return False
        try:
            if diff(expr, x) == nan or not Eq(diff(expr, x, x), 0):
                return False
        except TypeError:
            return False
    return True

def __init__(self, exprs):
        """
        Initialize a Scope, which represents a group of :class:`Access` objects
        extracted from some expressions ``exprs``. The expressions are to be
        provided as they appear in program order.
        """
        exprs = as_tuple(exprs)
        assert all(isinstance(i, Eq) for i in exprs)

        self.reads = {}
        self.writes = {}
        for i, e in enumerate(exprs):
            # reads
            for j in retrieve_indexed(e.rhs):
                v = self.reads.setdefault(j.base.function, [])
                mode = 'R' if not q_inc(e) else 'RI'
                v.append(TimedAccess(j, mode, i))
            # write
            if e.lhs.is_Indexed:
                v = self.writes.setdefault(e.lhs.base.function, [])
                mode = 'W' if not q_inc(e) else 'WI'
                v.append(TimedAccess(e.lhs, mode, i))

def __init__(self, expr, dtype=None):
        assert isinstance(expr, Eq)
        assert isinstance(expr.lhs, (Symbol, Indexed))
        self.expr = expr
        self.dtype = dtype

        # Traverse /expression/ to determine meta information
        # Note: at this point, expressions have already been indexified
        self.reads = [i for i in retrieve_terminals(self.expr.rhs)
                      if isinstance(i, (types.Indexed, types.Symbol))]
        self.reads = filter_ordered(self.reads)
        self.functions = [self.write] + [i.base.function for i in self.reads]
        self.functions = filter_ordered(self.functions)
        # Filter collected dimensions and functions
        self.dimensions = flatten(i.indices for i in self.functions)
        self.dimensions = filter_ordered(self.dimensions)

def execute_devito(ui, spacing=0.01, a=0.5, timesteps=500):
    """Execute diffusion stencil using the devito Operator API."""
    nx, ny = ui.shape
    dx2, dy2 = spacing**2, spacing**2
    dt = dx2 * dy2 / (2 * a * (dx2 + dy2))
    # Allocate the grid and set initial condition
    # Note: This should be made simpler through the use of defaults
    grid = Grid(shape=(nx, ny))
    u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
    u.data[0, :] = ui[:]

    # Derive the stencil according to devito conventions
    eqn = Eq(u.dt, a * (u.dx2 + u.dy2))
    stencil = sympy.solve(eqn, u.forward)[0]
    op = Operator(Eq(u.forward, stencil))

    # Execute the generated Devito stencil operator
    tstart = time.time()
    op.apply(u=u, t=timesteps, dt=dt)
    runtime = time.time() - tstart
    log("Devito: Diffusion with dx=%0.4f, dy=%0.4f, executed %d timesteps in %f seconds"
        % (spacing, spacing, timesteps, runtime))
    return u.data[1, :], runtime

def plot():
    e = Symbol('e')
    y = Symbol('y')
    n = Symbol('n')
    generalized_vc_bounds = (original_vc_bound, rademacher_penalty_bound)
    growth_function_bound = generate_growth_function_bound(50)
    p1 = plot(original_vc_bound(n, 0.05, growth_function_bound), (n,100, 15000), show=False, line_color = 'black')
    p2 = plot(rademacher_penalty_bound(n, 0.05, growth_function_bound), (n,100, 15000), show=False, line_color = 'blue')
    plot_implicit(Eq(e, parrondo_van_den_broek_right(e, n, 0.05, growth_function_bound)), (n,100, 15000), (e,0,5))
    # plot_implicit(Eq(e, devroye(e, n, 0.05, growth_function_bound)), (n,100, 1000), (e,0,5))
    p1.extend(p2)
    p1.show()

# BIAS AND VARIANCE

def test_density():
    x = Symbol('x')
    l = Symbol('l', positive=True)
    rate = Beta(l, 2, 3)
    X = Poisson(x, rate)
    assert isinstance(pspace(X), ProductPSpace)
    assert density(X, Eq(rate, rate.symbol)) == PoissonDistribution(l)

def test_signsimp():
    e = x*(-x + 1) + x*(x - 1)
    assert signsimp(Eq(e, 0)) is S.true
    assert Abs(x - 1) == Abs(1 - x)

def runtest_autowrap_matrix_vector(language, backend):
    has_module('numpy')
    x, y = symbols('x y', cls=IndexedBase)
    expr = Eq(y[i], A[i, j]*x[j])
    mv = autowrap(expr, language, backend)

    # compare with numpy's dot product
    M = numpy.random.rand(10, 20)
    x = numpy.random.rand(20)
    y = numpy.dot(M, x)
    assert numpy.sum(numpy.abs(y - mv(M, x))) < 1e-13

def runtest_autowrap_matrix_matrix(language, backend):
    has_module('numpy')
    expr = Eq(C[i, j], A[i, k]*B[k, j])
    matmat = autowrap(expr, language, backend)

    # compare with numpy's dot product
    M1 = numpy.random.rand(10, 20)
    M2 = numpy.random.rand(20, 15)
    M3 = numpy.dot(M1, M2)
    assert numpy.sum(numpy.abs(M3 - matmat(M1, M2))) < 1e-13

def test_dict_from_expr():
    dict_from_expr(Eq(x, 1)) == ({(1,): Integer(1)}, (x,))
    raises(PolynomialError, lambda: dict_from_expr(A*B - B*A))

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 main():
    print("Hydrogen radial wavefunctions:")
    a, r = symbols("a r")
    print("R_{21}:")
    pprint(R_nl(2, 1, a, r))
    print("R_{60}:")
    pprint(R_nl(6, 0, a, r))

    print("Normalization:")
    i = Integral(R_nl(1, 0, 1, r)**2 * r**2, (r, 0, oo))
    pprint(Eq(i, i.doit()))
    i = Integral(R_nl(2, 0, 1, r)**2 * r**2, (r, 0, oo))
    pprint(Eq(i, i.doit()))
    i = Integral(R_nl(2, 1, 1, r)**2 * r**2, (r, 0, oo))
    pprint(Eq(i, i.doit()))

def p_booleqExpr_impl(self,p):
        '''booleqExpr : relationExpr BOOLEQ relationExpr
                      | relationExpr NEQ    relationExpr
        '''
        if p[2] == "NEQ":
            boolOp = sympy.Ne
        else:
            boolOp = sympy.Eq
        self.checkExactUnits(p[1].astUnit,p[3].astUnit)
        p[0] = AST(boolOp(p[1].sympy,p[3].sympy),self.boolean())

def test_matrix_to_expr_operator(matrix_form, expr1, expr2):
    assert sp.simplify(
        sp.Eq(matrix_to_expr_operator(matrix_form)(expr1), expr2)
    )

def test_matrix_to_expr_operator_double_eval(matrix_form, expr1, expr2):
    expr_operator = matrix_to_expr_operator(matrix_form)
    assert sp.simplify(sp.Eq(expr_operator(expr1), expr2))
    assert sp.simplify(sp.Eq(expr_operator(expr1), expr2))

def test_probability_rewrite():
    X = Normal('X', 2, 3)
    Y = Normal('Y', 3, 4)
    Z = Poisson('Z', 4)
    W = Poisson('W', 3)
    x, y, w, z = symbols('x, y, w, z')

    assert Variance(w).rewrite(Expectation) == 0
    assert Variance(X).rewrite(Expectation) == Expectation(X ** 2) - Expectation(X) ** 2
    assert Variance(X, condition=Y).rewrite(Expectation) == Expectation(X ** 2, Y) - Expectation(X, Y) ** 2
    assert Variance(X, Y) != Expectation(X**2) - Expectation(X)**2
    assert Variance(X + z).rewrite(Expectation) == Expectation((X + z) ** 2) - Expectation(X + z) ** 2
    assert Variance(X * Y).rewrite(Expectation) == Expectation(X ** 2 * Y ** 2) - Expectation(X * Y) ** 2

    assert Covariance(w, X).rewrite(Expectation) == -w*Expectation(X) + Expectation(w*X)
    assert Covariance(X, Y).rewrite(Expectation) == Expectation(X*Y) - Expectation(X)*Expectation(Y)
    assert Covariance(X, Y, condition=W).rewrite(Expectation) == Expectation(X * Y, W) - Expectation(X, W) * Expectation(Y, W)

    w, x, z = symbols("W, x, z")
    px = Probability(Eq(X, x))
    pz = Probability(Eq(Z, z))

    assert Expectation(X).rewrite(Probability) == Integral(x*px, (x, -oo, oo))
    assert Expectation(Z).rewrite(Probability) == Sum(z*pz, (z, 0, oo))
    assert Variance(X).rewrite(Probability) == Integral(x**2*px, (x, -oo, oo)) - Integral(x*px, (x, -oo, oo))**2
    assert Variance(Z).rewrite(Probability) == Sum(z**2*pz, (z, 0, oo)) - Sum(z*pz, (z, 0, oo))**2

    assert Variance(X, condition=Y).rewrite(Probability) == Integral(x**2*Probability(Eq(X, x), Y), (x, -oo, oo)) - \
                                                            Integral(x*Probability(Eq(X, x), Y), (x, -oo, oo))**2

def test_density():
    x = Symbol('x')
    l = Symbol('l', positive=True)
    rate = Beta(l, 2, 3)
    X = Poisson(x, rate)
    assert isinstance(pspace(X), ProductPSpace)
    assert density(X, Eq(rate, rate.symbol)) == PoissonDistribution(l)

def test_euler_interface():
    x = Function('x')
    y = Symbol('y')
    t = Symbol('t')
    raises(TypeError, lambda: euler())
    raises(TypeError, lambda: euler(D(x(t), t)*y(t), [x(t), y]))
    raises(ValueError, lambda: euler(D(x(t), t)*x(y), [x(t), x(y)]))
    raises(TypeError, lambda: euler(D(x(t), t)**2, x(0)))
    assert euler(D(x(t), t)**2/2, {x(t)}) == [Eq(-D(x(t), t, t))]
    assert euler(D(x(t), t)**2/2, x(t), {t}) == [Eq(-D(x(t), t, t))]