def chapman_vec(z_vec, Nm_vec, Hm_vec, H_O_vec):
"""
Vectorized implementation of the Chapman function evaluation
routine :func:`chapman`. The input arguments must be sequences
with the same length and the output is an :class:`NP.ndarray` with
that length.
"""
try:
chapman_vec._chapman_sym_f
except AttributeError:
sym_vars = SYM.symbols('z Nm Hm H_O')
chapman_vec._chapman_sym_f = SYM.lambdify(sym_vars,
chapman_sym(*sym_vars),
modules='numexpr')
return chapman_vec._chapman_sym_f(z_vec,
Nm_vec,
Hm_vec,
H_O_vec)
python类lambdify()的实例源码
def test_clambdify():
x = Symbol('x')
y = Symbol('y')
z = Symbol('z')
f1 = sqrt(x*y)
pf1 = lambdify((x, y), f1, 'math')
cf1 = clambdify((x, y), f1)
for i in xrange(10):
if cf1(i, 10 - i) != pf1(i, 10 - i):
raise ValueError("Values should be equal")
f2 = (x - y) / z * pi
pf2 = lambdify((x, y, z), f2, 'math')
cf2 = clambdify((x, y, z), f2)
if round(pf2(1, 2, 3), 14) != round(cf2(1, 2, 3), 14):
raise ValueError("Values should be equal")
# FIXME: slight difference in precision
def __init__(self, settings):
Trajectory.__init__(self, settings)
# setup symbolic expressions
tau, k = sp.symbols('tau, k')
gamma = self._settings["differential_order"] + 1
alpha = sp.factorial(2 * gamma + 1)
f = sp.binomial(gamma, k) * (-1) ** k * tau ** (gamma + k + 1) / (gamma + k + 1)
phi = alpha / sp.factorial(gamma) ** 2 * sp.summation(f, (k, 0, gamma))
# differentiate phi(tau), index in list corresponds to order
dphi_sym = [phi] # init with phi(tau)
for order in range(self._settings["differential_order"]):
dphi_sym.append(dphi_sym[-1].diff(tau))
# lambdify
self.dphi_num = []
for der in dphi_sym:
self.dphi_num.append(sp.lambdify(tau, der, 'numpy'))
def __init__(self, *args, **kwargs):
self.args = args
self.f, self.gradient = None, ColorGradient()
if len(args) == 1 and not isinstance(args[0], Basic) and callable(args[0]):
self.f = args[0]
elif len(args) == 1 and isinstance(args[0], str):
if args[0] in default_color_schemes:
cs = default_color_schemes[args[0]]
self.f, self.gradient = cs.f, cs.gradient.copy()
else:
self.f = lambdify('x,y,z,u,v', args[0])
else:
self.f, self.gradient = self._interpret_args(args, kwargs)
self._test_color_function()
if not isinstance(self.gradient, ColorGradient):
raise ValueError("Color gradient not properly initialized. "
"(Not a ColorGradient instance.)")
def test_imps():
# Here we check if the default returned functions are anonymous - in
# the sense that we can have more than one function with the same name
f = implemented_function('f', lambda x: 2*x)
g = implemented_function('f', lambda x: math.sqrt(x))
l1 = lambdify(x, f(x))
l2 = lambdify(x, g(x))
assert str(f(x)) == str(g(x))
assert l1(3) == 6
assert l2(3) == math.sqrt(3)
# check that we can pass in a Function as input
func = sympy.Function('myfunc')
assert not hasattr(func, '_imp_')
my_f = implemented_function(func, lambda x: 2*x)
assert hasattr(func, '_imp_')
# Error for functions with same name and different implementation
f2 = implemented_function("f", lambda x: x + 101)
raises(ValueError, lambda: lambdify(x, f(f2(x))))
def test_dummification():
t = symbols('t')
F = Function('F')
G = Function('G')
#"\alpha" is not a valid python variable name
#lambdify should sub in a dummy for it, and return
#without a syntax error
alpha = symbols(r'\alpha')
some_expr = 2 * F(t)**2 / G(t)
lam = lambdify((F(t), G(t)), some_expr)
assert lam(3, 9) == 2
lam = lambdify(sin(t), 2 * sin(t)**2)
assert lam(F(t)) == 2 * F(t)**2
#Test that \alpha was properly dummified
lam = lambdify((alpha, t), 2*alpha + t)
assert lam(2, 1) == 5
raises(SyntaxError, lambda: lambdify(F(t) * G(t), F(t) * G(t) + 5))
raises(SyntaxError, lambda: lambdify(2 * F(t), 2 * F(t) + 5))
raises(SyntaxError, lambda: lambdify(2 * F(t), 4 * F(t) + 5))
def test_special_printers():
class IntervalPrinter(LambdaPrinter):
"""Use ``lambda`` printer but print numbers as ``mpi`` intervals. """
def _print_Integer(self, expr):
return "mpi('%s')" % super(IntervalPrinter, self)._print_Integer(expr)
def _print_Rational(self, expr):
return "mpi('%s')" % super(IntervalPrinter, self)._print_Rational(expr)
def intervalrepr(expr):
return IntervalPrinter().doprint(expr)
expr = sympy.sqrt(sympy.sqrt(2) + sympy.sqrt(3)) + sympy.S(1)/2
func0 = lambdify((), expr, modules="mpmath", printer=intervalrepr)
func1 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter)
func2 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter())
mpi = type(mpmath.mpi(1, 2))
assert isinstance(func0(), mpi)
assert isinstance(func1(), mpi)
assert isinstance(func2(), mpi)
def application():
"""makes list of the difference between the estimated derivative for h=0.01 and
the analytical derivative"""
application_list = []
x = symbols('x')
print
print '%20s %20s' % ('Function', 'Error')
print
funcvar_tuples = [(lambda t: exp(t), 'exp(x)', float(0)), (lambda t : exp(-2*t**2), \
'exp(-2*x**2)', float(0)), (lambda t: cos(t), 'cos(x)', 2*pi), (lambda t: log(t), \
'log(x)', float(1))]
for i in funcvar_tuples:
prime = diff(i[1], x)
exact = lambdify([x], prime)
error = numdiff(i[0], i[2], 0.01) - exact(i[2])
print '%20s %21.2E' % (i[1], error)
print
application_list.append(error)
return application_list
def __init__(self, *args, **kwargs):
self.args = args
self.f, self.gradient = None, ColorGradient()
if len(args) == 1 and not isinstance(args[0], Basic) and callable(args[0]):
self.f = args[0]
elif len(args) == 1 and isinstance(args[0], str):
if args[0] in default_color_schemes:
cs = default_color_schemes[args[0]]
self.f, self.gradient = cs.f, cs.gradient.copy()
else:
self.f = lambdify('x,y,z,u,v', args[0])
else:
self.f, self.gradient = self._interpret_args(args, kwargs)
self._test_color_function()
if not isinstance(self.gradient, ColorGradient):
raise ValueError("Color gradient not properly initialized. "
"(Not a ColorGradient instance.)")
def test_sum():
if not np:
skip("NumPy not installed")
s = Sum(x ** i, (i, a, b))
f = lambdify((a, b, x), s, 'numpy')
a_, b_ = 0, 10
x_ = np.linspace(-1, +1, 10)
assert np.allclose(f(a_, b_, x_), sum(x_ ** i_ for i_ in range(a_, b_ + 1)))
s = Sum(i * x, (i, a, b))
f = lambdify((a, b, x), s, 'numpy')
a_, b_ = 0, 10
x_ = np.linspace(-1, +1, 10)
assert np.allclose(f(a_, b_, x_), sum(i_ * x_ for i_ in range(a_, b_ + 1)))
def test_mod():
if not np:
skip("NumPy not installed")
e = Mod(a, b)
f = lambdify((a, b), e)
a_ = np.array([0, 1, 2, 3])
b_ = 2
assert np.array_equal(f(a_, b_), [0, 1, 0, 1])
a_ = np.array([0, 1, 2, 3])
b_ = np.array([2, 2, 2, 2])
assert np.array_equal(f(a_, b_), [0, 1, 0, 1])
a_ = np.array([2, 3, 4, 5])
b_ = np.array([2, 3, 4, 5])
assert np.array_equal(f(a_, b_), [0, 0, 0, 0])
def test_sum():
# In each case, test eval() the lambdarepr() to make sure that
# it evaluates to the same results as the symbolic expression
s = Sum(x ** i, (i, a, b))
l = lambdarepr(s)
assert l == "(builtins.sum(x**i for i in range(a, b+1)))"
assert (lambdify((x, a, b), s)(2, 3, 8) ==
s.subs([(x, 2), (a, 3), (b, 8)]).doit())
s = Sum(i * x, (i, a, b))
l = lambdarepr(s)
assert l == "(builtins.sum(i*x for i in range(a, b+1)))"
assert (lambdify((x, a, b), s)(2, 3, 8) ==
s.subs([(x, 2), (a, 3), (b, 8)]).doit())
def __missing__(self, i):
args = ['p%d'%d for d in range(i+1)]
f = green(self._symfunc, BezierCurve[i])
f = sp.gcd_terms(f.collect(sum(P,()))) # Optimize
return sp.lambdify(args, f)
def __init__(self, stn_pos, **kwds):
"""
Setup Chapman electron density profile slant integrator.
"""
z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
f_sym = chapman_sym_scaled(z, Nm, Hm, H_O)
f = SYM.lambdify((z, Nm, Hm, H_O),
f_sym,
modules='numexpr')
wrapper = lambda pos, *args: f(pos.height / 1e3, *args)
super(ChapmanSI, self).__init__(wrapper, stn_pos, **kwds)
def __init__(self, stn_pos, **kwds):
"""
Setup Chapman electron density profile F2 peak density derivative
slant integrator.
"""
z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
f_sym = chapman_sym_scaled(z, Nm, Hm, H_O)
DNm_sym = SYM.diff(f_sym, Nm)
f = SYM.lambdify((z, Hm, H_O),
DNm_sym,
modules='numexpr')
wrapper = lambda pos, *args: f(pos.height / 1e3, *args)
super(DNmChapmanSI, self).__init__(wrapper, stn_pos, **kwds)
def __init__(self, stn_pos, **kwds):
"""
Setup Chapman electron density profile F2 peak height derivative
slant integrator.
"""
z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
f_sym = chapman_sym_scaled(z, Nm, Hm, H_O)
DHm_sym = SYM.diff(f_sym, Hm)
f = SYM.lambdify((z, Nm, Hm, H_O),
DHm_sym,
modules='numexpr')
wrapper = lambda pos, *args: f(pos.height / 1e3, *args)
super(DHmChapmanSI, self).__init__(wrapper, stn_pos, **kwds)
variable_stepsize_custom_constraints.py 文件源码
项目:toppra
作者: hungpham2511
项目源码
文件源码
阅读 19
收藏 0
点赞 0
评论 0
def evaluate_constraint(expressions, ss):
"""Evaluate the canonical constraints `expressions` over the
discretized interval `ss`.
Parameters
----------
Returns
-------
"""
m = len(expressions)
N = len(ss) - 1
a = np.zeros((N + 1, m))
b = np.zeros((N + 1, m))
c = np.zeros((N + 1, m))
for index, e in enumerate(expressions):
if e.rel_op == "<=":
canon_form = e.lhs - e.rhs
elif e.rel_op == ">=":
canon_form = e.rhs - e.lhs
else:
raise ValueError("Relation Operation need to be `<=` or `>=`.")
f = sympy.lambdify(sympy.symbols('u, x, s'), canon_form)
for i in xrange(N + 1):
a[i, index] = f(1, 0, ss[i]) - f(0, 0, ss[i])
b[i, index] = f(0, 1, ss[i]) - f(0, 0, ss[i])
c[i, index] = f(0, 0, ss[i])
return a, b, c
# Define symbolic expression and evaluation
def vectors(self, mesh):
""" Get Vectors sig, sr. jx from sympy"""
j, Sig = self.fcts()
f_jr = sympy.lambdify((r, z), j[0], 'numpy')
f_jz = sympy.lambdify((r, z), j[1], 'numpy')
f_sigr = sympy.lambdify((r, z), Sig[0], 'numpy')
# f_sigz = sympy.lambdify((r,z), Sig[1], 'numpy')
jr = f_jr(mesh.gridFx[:, 0], mesh.gridFx[:, 2])
jz = f_jz(mesh.gridFz[:, 0], mesh.gridFz[:, 2])
sigr = f_sigr(mesh.gridCC[:, 0], mesh.gridCC[:, 2])
return sigr, np.r_[jr, jz]
def vectors(self, mesh):
""" Get Vectors sig, sr. jx from sympy"""
j, Sig = self.fcts()
f_jr = sympy.lambdify((r, z), j[0], 'numpy')
f_jz = sympy.lambdify((r, z), j[1], 'numpy')
f_sigr = sympy.lambdify((r, z), Sig[0], 'numpy')
f_sigz = sympy.lambdify((r, z), Sig[3], 'numpy')
jr = f_jr(mesh.gridFx[:, 0], mesh.gridFx[:, 2])
jz = f_jz(mesh.gridFz[:, 0], mesh.gridFz[:, 2])
sigr = f_sigr(mesh.gridCC[:, 0], mesh.gridCC[:, 2])
sigz = f_sigz(mesh.gridCC[:, 0], mesh.gridCC[:, 2])
return np.c_[sigr, sigr, sigz], np.r_[jr, jz]
def vectors(self, mesh):
""" Get Vectors sig, sr. jx from sympy"""
h, Sig = self.fcts()
f_h = sympy.lambdify((r, z), h[0], 'numpy')
f_sig = sympy.lambdify((r, z), Sig[0], 'numpy')
ht = f_h(mesh.gridEy[:, 0], mesh.gridEy[:, 2])
sig = f_sig(mesh.gridCC[:, 0], mesh.gridCC[:, 2])
return sig, np.r_[ht]
def vectors(self, mesh):
h, Sig = self.fcts()
f_h = sympy.lambdify((r, z), h[0], 'numpy')
f_sig = sympy.lambdify((r, z), Sig[0], 'numpy')
ht = f_h(mesh.gridEy[:, 0], mesh.gridEy[:, 2])
sig = f_sig(mesh.gridCC[:, 0], mesh.gridCC[:, 2])
return np.c_[sig, sig, sig], np.r_[ht]
def test_tanh_sinh(f, a, b, exact):
# test fine error estimate
mp.dps = 50
tol = 10**(-mp.dps)
tol2 = 10**(-mp.dps+1)
t = sympy.Symbol('t')
f_derivatives = {
1: sympy.lambdify(t, sympy.diff(f(t), t, 1), modules=['mpmath']),
2: sympy.lambdify(t, sympy.diff(f(t), t, 2), modules=['mpmath']),
}
value, _ = quadpy.line_segment.tanh_sinh(
f, a, b, tol,
f_derivatives=f_derivatives
)
assert abs(value - exact) < tol2
# test with crude estimate
value, _ = quadpy.line_segment.tanh_sinh(f, a, b, tol)
assert abs(value - exact) < tol2
return
# Test functions with singularities at both ends.
def test_singularities_at_both_ends(f_left, f_right, b, exact):
# test fine error estimate
tol = 10**(-mp.dps)
t = sympy.Symbol('t')
fl = {
0: f_left,
1: sympy.lambdify(t, sympy.diff(f_left(t), t, 1), modules=['mpmath']),
2: sympy.lambdify(t, sympy.diff(f_left(t), t, 2), modules=['mpmath']),
}
fr = {
0: f_right,
1: sympy.lambdify(t, sympy.diff(f_right(t), t, 1), modules=['mpmath']),
2: sympy.lambdify(t, sympy.diff(f_right(t), t, 2), modules=['mpmath']),
}
value, _ = quadpy.line_segment.tanh_sinh_lr(fl, fr, b, tol)
tol2 = 10**(-mp.dps+1)
assert abs(value - exact) < tol2
# # test with crude estimate
# fl = {0: f_left}
# fr = {0: f_right}
# value, _ = quadpy.line_segment.tanh_sinh_lr(fl, fr, b, tol)
# tol2 = 10**(-mp.dps+2)
# assert abs(value - exact) < tol2
return
def __init__(self, settings):
settings.update(output_dim=4)
super().__init__(settings)
params = sp.symbols('x1, x2, x3, x4, tau')
x1, x2, x3, x4, tau = params
x = [x1, x2, x3, x4]
h = sp.Matrix([[x1]])
f = sp.Matrix([[x2],
[st.B * x1 * x4 ** 2 - st.B * st.G * sin(x3)],
[x4],
[(tau - 2 * st.M * x1 * x2 * x4
- st.M * st.G * x1 * cos(x3)) / (
st.M * x1 ** 2 + st.J + st.Jb)]])
q = sp.Matrix(pm.lie_derivatives(h, f, x, len(x) - 1))
dq = q.jacobian(x)
if dq.rank() != len(x):
raise Exception("System might not be observable")
# gets p = [p0, p1, ... pn-1]
p = pm.char_coefficients(self._settings["poles"])
k = p[::-1]
l = dq.inv() @ k
mat2array = [{'ImmutableMatrix': np.array}, 'numpy']
self.h_func = sp.lambdify((x1, x2, x3, x4, tau), h, modules=mat2array)
self.l_func = sp.lambdify((x1, x2, x3, x4, tau), l, modules=mat2array)
self.f_func = sp.lambdify((x1, x2, x3, x4, tau), f, modules=mat2array)
self.output = np.array(self._settings["initial state"], dtype=float)
def __init__(self, settings):
Trajectory.__init__(self, settings)
# calculate symbolic derivatives up to order n
t, a, f, off, p = sp.symbols("t, a, f, off, p")
self.yd_sym = []
harmonic = a * sp.sin(2 * sp.pi * f * t + p) + off
for idx in range(settings["differential_order"] + 1):
self.yd_sym.append(harmonic.diff(t, idx))
# lambdify
for idx, val in enumerate(self.yd_sym):
self.yd_sym[idx] = sp.lambdify((t, a, f, off, p), val, "numpy")
def lambdify_expr(expr):
return None_on_RuntimeError(lambdify(t, expr, scipy_translations_autoeye, printer=MatrixNumPyPrinter))
def mean_log10_relative_error(exact, approx):
import numpy as np
return np.mean(np.log10(abs(relative_error(exact, approx))))
# Function to create plot like in "Computing the Matrix Exponential in Burnup
# Calculations", Pusa and Leppa?nen:
# mpmath.cplot(lambdify(t, rat_func14 - exp(-t), 'mpmath'), re=[0, 100], im=[-30, 30], color=lambda i: -mpmath.floor(mpmath.log(abs(i), 10))/(30 - mpmath.floor(mpmath.log(abs(i), 10))), points=100000, verbose=True)
def cplot_in_terminal(expr, *args, prec=None, logname=None, color=lambda i:
-mpmath.floor(mpmath.log(abs(i), 10))/(30 -
mpmath.floor(mpmath.log(abs(i), 10))), points=1000000, **kwargs):
"""
Run mpmath.cplot() but show in terminal if possible
"""
kwargs['color'] = color
kwargs['points'] = points
from mpmath import cplot
if prec:
mpmath.mp.dps = prec
f = lambdify(t, expr, mpmath)
try:
from iterm2_tools.images import display_image_bytes
except ImportError:
if logname:
os.makedirs('plots', exist_ok=True)
file = 'plots/%s.png' % logname
else:
file = None
cplot(f, *args, file=file, **kwargs)
else:
from io import BytesIO
b = BytesIO()
cplot(f, *args, **kwargs, file=b)
if logname:
os.makedirs('plots', exist_ok=True)
with open('plots/%s.png' % logname, 'wb') as f:
f.write(b.getvalue())
print(display_image_bytes(b.getvalue()))
def plot_in_terminal(expr, *args, prec=None, logname=None, file=None, **kwargs):
"""
Run mpmath.plot() but show in terminal if possible
"""
from mpmath import plot
if logname:
os.makedirs('plots', exist_ok=True)
file = 'plots/%s.png' % logname
if prec:
mpmath.mp.dps = prec
if isinstance(expr, (list, tuple)):
f = [lambdify(t, i, 'mpmath') for i in expr]
else:
f = lambdify(t, expr, 'mpmath')
try:
if hasattr(sys.stdout, 'isatty') and sys.stdout.isatty():
from iterm2_tools.images import display_image_bytes
else:
raise ImportError
except ImportError:
plot(f, *args, file=file, **kwargs)
else:
# mpmath.plot ignores the axes argument if file is given, so let
# file=False, disable this.
if 'axes' in kwargs:
file=False
if file is not False:
from io import BytesIO
b = BytesIO()
else:
b = None
plot(f, *args, **kwargs, file=b)
if file:
with open(file, 'wb') as f:
f.write(b.getvalue())
if b:
print(display_image_bytes(b.getvalue()))
def _compile(expr, variable=mu):
"""compiles a sympy expression"""
return lambdify(variable, expr) if isinstance(expr, Basic) else expr