我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用sympy.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_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))]
def test_euler_pendulum(): x = Function('x') t = Symbol('t') L = D(x(t), t)**2/2 + cos(x(t)) assert euler(L, x(t), t) == [Eq(-sin(x(t)) - D(x(t), t, t))]
def test_euler_sineg(): psi = Function('psi') t = Symbol('t') x = Symbol('x') L = D(psi(t, x), t)**2/2 - D(psi(t, x), x)**2/2 + cos(psi(t, x)) assert euler(L, psi(t, x), [t, x]) == [Eq(-sin(psi(t, x)) - D(psi(t, x), t, t) + D(psi(t, x), x, x))]
def test_preview_latex_construct_in_expr(): # see PR 9801 x = Symbol('x') pw = Piecewise((1, Eq(x, 0)), (0, True)) obj = BytesIO() try: preview(pw, output='png', viewer='BytesIO', outputbuffer=obj) except RuntimeError: pass # latex not installed on CI server
def test_ccode_Relational(): from sympy import Eq, Ne, Le, Lt, Gt, Ge assert ccode(Eq(x, y)) == "x == y" assert ccode(Ne(x, y)) == "x != y" assert ccode(Le(x, y)) == "x <= y" assert ccode(Lt(x, y)) == "x < y" assert ccode(Gt(x, y)) == "x > y" assert ccode(Ge(x, y)) == "x >= y"
def test_rcode_Relational(): from sympy import Eq, Ne, Le, Lt, Gt, Ge assert rcode(Eq(x, y)) == "x == y" assert rcode(Ne(x, y)) == "x != y" assert rcode(Le(x, y)) == "x <= y" assert rcode(Lt(x, y)) == "x < y" assert rcode(Gt(x, y)) == "x > y" assert rcode(Ge(x, y)) == "x >= y"
def test_rcode_Indexed_without_looking_for_contraction(): len_y = 5 y = IndexedBase('y', shape=(len_y,)) x = IndexedBase('x', shape=(len_y,)) Dy = IndexedBase('Dy', shape=(len_y-1,)) i = Idx('i', len_y-1) e=Eq(Dy[i], (y[i+1]-y[i])/(x[i+1]-x[i])) code0 = rcode(e.rhs, assign_to=e.lhs, contract=False) assert code0 == 'Dy[i] = (y[%s] - y[i])/(x[%s] - x[i]);' % (i + 1, i + 1)
def solveEquation(expr): expr = expr.replace('^', '**') members = expr.split('=') if len(members) != 2: raise BadFormatException('Bad number of equals.') from sympy.abc import * eq = sympy.Eq(*map(eval, members)) return [{repr(j): repr(k) for j, k in i.items()} for i in _ensureList(sympy.solve(eq, dict=1))]
def test_one_variable(self): """ It can eliminate ``x`` from the system ``(a = x, b = x)``. """ equations = [sympy.Eq(a, x), sympy.Eq(b, x)] new_equations = vadouvan.symbolic_utils.eliminate(equations, [x]) assert len(new_equations) == 1 assert new_equations[0].free_symbols == {a, b}
def test_two_variables(self): """ It can eliminate ``x, y`` from the system ``(a = 3x - 4y, b = 4x + 5y, c = 6x + 2y)``. """ equations = [ sympy.Eq(a, 3 * x - 4 * y), sympy.Eq(b, 4 * x + 5 * y), sympy.Eq(c, 6 * x + 2 * y), ] new_equations = vadouvan.symbolic_utils.eliminate(equations, [x, y]) assert len(new_equations) == 1 assert new_equations[0].free_symbols == {a, b, c}
def test_many_variables(self): """ It can eliminate all x's from the chain ``(a = x0, x0 = x1, ..., x19 = x20, x20 = b)``. """ symbols = list(itertools.islice(sympy.numbered_symbols("x"), 20)) equations = [sympy.Eq(symbols[i], symbols[i + 1]) for i in range(len(symbols) - 1)] equations.append(sympy.Eq(a, symbols[0])) equations.append(sympy.Eq(b, symbols[-1])) new_equations = vadouvan.symbolic_utils.eliminate(equations, symbols) assert len(new_equations) == 1 assert new_equations[0].free_symbols == {a, b}
def test_underconstrained_2(self): """ Underconstrained system returns an empty collection. """ equations = [sympy.Eq(x, y)] new_equations = vadouvan.symbolic_utils.eliminate(equations, [x, y]) assert len(new_equations) == 0
def odes(self): """Return a list of differential equations describing the evolution of the species concentrations in time.""" t = sp.Symbol('t', real = True) x = [sp.Function(s) for s in self.species] derivs = self.equations() return [sp.Eq(sp.Derivative(x[j](t), t), derivs[j].subs([(sp.Symbol(self.species[i]), x[i](t)) for i in range(self.n_species)])) for j in range(self.n_species)]
def __new__(cls, *args, **kwargs): kwargs['evaluate'] = False obj = sympy.Eq.__new__(cls, *args, **kwargs) return obj
def ccode(expr, **settings): """Generate C++ code from an expression calling CodePrinter class :param expr: The expression :param settings: A dictionary of settings for code printing :returns: The resulting code as a string. If it fails, then it returns the expr """ if isinstance(expr, Eq): return ccode_eq(expr) try: return CodePrinter(settings).doprint(expr, None) except: return expr
def first_touch(array): """ Uses an Operator to initialize the given array in the same pattern that would later be used to access it. """ devito.Operator(Eq(array, 0.))()
def execute_lambdify(ui, spacing=0.01, a=0.5, timesteps=500): """Execute diffusion stencil using vectorised numpy array accesses.""" nx, ny = ui.shape dx2, dy2 = spacing**2, spacing**2 dt = dx2 * dy2 / (2 * a * (dx2 + dy2)) u = np.concatenate((ui, np.zeros_like(ui))).reshape((2, nx, ny)) def diffusion_stencil(): """Create stencil and substitutions for the diffusion equation""" p = sympy.Function('p') x, y, t, h, s = sympy.symbols('x y t h s') dx2 = p(x, y, t).diff(x, x).as_finite_difference([x - h, x, x + h]) dy2 = p(x, y, t).diff(y, y).as_finite_difference([y - h, y, y + h]) dt = p(x, y, t).diff(t).as_finite_difference([t, t + s]) eqn = sympy.Eq(dt, a * (dx2 + dy2)) stencil = sympy.solve(eqn, p(x, y, t + s))[0] return stencil, (p(x, y, t), p(x + h, y, t), p(x - h, y, t), p(x, y + h, t), p(x, y - h, t), s, h) stencil, subs = diffusion_stencil() kernel = sympy.lambdify(subs, stencil, 'numpy') # Execute timestepping loop with alternating buffers tstart = time.time() for ti in range(timesteps): t0 = ti % 2 t1 = (ti + 1) % 2 u[t1, 1:-1, 1:-1] = kernel(u[t0, 1:-1, 1:-1], u[t0, 2:, 1:-1], u[t0, :-2, 1:-1], u[t0, 1:-1, 2:], u[t0, 1:-1, :-2], dt, spacing) runtime = time.time() - tstart log("Lambdify: Diffusion with dx=%0.4f, dy=%0.4f, executed %d timesteps in %f seconds" % (spacing, spacing, timesteps, runtime)) return u[ti % 2, :, :], runtime
def test_unpolarify(): from sympy import (exp_polar, polar_lift, exp, unpolarify, sin, principal_branch) from sympy import gamma, erf, sin, tanh, uppergamma, Eq, Ne from sympy.abc import x p = exp_polar(7*I) + 1 u = exp(7*I) + 1 assert unpolarify(1) == 1 assert unpolarify(p) == u assert unpolarify(p**2) == u**2 assert unpolarify(p**x) == p**x assert unpolarify(p*x) == u*x assert unpolarify(p + x) == u + x assert unpolarify(sqrt(sin(p))) == sqrt(sin(u)) # Test reduction to principal branch 2*pi. t = principal_branch(x, 2*pi) assert unpolarify(t) == x assert unpolarify(sqrt(t)) == sqrt(t) # Test exponents_only. assert unpolarify(p**p, exponents_only=True) == p**u assert unpolarify(uppergamma(x, p**p)) == uppergamma(x, p**u) # Test functions. assert unpolarify(sin(p)) == sin(u) assert unpolarify(tanh(p)) == tanh(u) assert unpolarify(gamma(p)) == gamma(u) assert unpolarify(erf(p)) == erf(u) assert unpolarify(uppergamma(x, p)) == uppergamma(x, p) assert unpolarify(uppergamma(sin(p), sin(p + exp_polar(0)))) == \ uppergamma(sin(u), sin(u + 1)) assert unpolarify(uppergamma(polar_lift(0), 2*exp_polar(0))) == \ uppergamma(0, 2) assert unpolarify(Eq(p, 0)) == Eq(u, 0) assert unpolarify(Ne(p, 0)) == Ne(u, 0) assert unpolarify(polar_lift(x) > 0) == (x > 0) # Test bools assert unpolarify(True) is True
def substitution_rule(integral): integrand, symbol = integral u_var = sympy.Dummy("u") substitutions = find_substitutions(integrand, symbol, u_var) if substitutions: ways = [] for u_func, c, substituted in substitutions: subrule = integral_steps(substituted, u_var) if contains_dont_know(subrule): continue if sympy.simplify(c - 1) != 0: _, denom = c.as_numer_denom() subrule = ConstantTimesRule(c, substituted, subrule, substituted, u_var) if denom.free_symbols: piecewise = [] could_be_zero = [] if isinstance(denom, sympy.Mul): could_be_zero = denom.args else: could_be_zero.append(denom) for expr in could_be_zero: if not sympy.ask(~sympy.Q.zero(expr)): substep = integral_steps(integrand.subs(expr, 0), symbol) if substep: piecewise.append(( substep, sympy.Eq(expr, 0) )) piecewise.append((subrule, True)) subrule = PiecewiseRule(piecewise, substituted, symbol) ways.append(URule(u_var, u_func, c, subrule, integrand, symbol)) if len(ways) > 1: return AlternativeRule(ways, integrand, symbol) elif ways: return ways[0] elif integrand.has(sympy.exp): u_func = sympy.exp(symbol) c = 1 substituted = integrand / u_func.diff(symbol) substituted = substituted.subs(u_func, u_var) if symbol not in substituted.free_symbols: return URule(u_var, u_func, c, integral_steps(substituted, u_var), integrand, symbol)
