Python sympy.lambdify函数代码示例

def test_numpy_matrix():
    if not numpy:
        skip("numpy not installed.")
    A = Matrix([[x, x*y], [sin(z) + 4, x**z]])
    sol_arr = numpy.array([[1, 2], [numpy.sin(3) + 4, 1]])
    #Lambdify array first, to ensure return to array as default
    f = lambdify((x, y, z), A, ['numpy'])
    numpy.testing.assert_allclose(f(1, 2, 3), sol_arr)
    #Check that the types are arrays and matrices
    assert isinstance(f(1, 2, 3), numpy.ndarray)

    # gh-15071
    class dot(Function):
        pass
    x_dot_mtx = dot(x, Matrix([[2], [1], [0]]))
    f_dot1 = lambdify(x, x_dot_mtx)
    inp = numpy.zeros((17, 3))
    assert numpy.all(f_dot1(inp) == 0)

    strict_kw = dict(allow_unknown_functions=False, inline=True, fully_qualified_modules=False)
    p2 = NumPyPrinter(dict(user_functions={'dot': 'dot'}, **strict_kw))
    f_dot2 = lambdify(x, x_dot_mtx, printer=p2)
    assert numpy.all(f_dot2(inp) == 0)

    p3 = NumPyPrinter(strict_kw)
    # The line below should probably fail upon construction (before calling with "(inp)"):
    raises(Exception, lambda: lambdify(x, x_dot_mtx, printer=p3)(inp))

def run_Lagrange_interp_abs_Cheb(N, ymin=None, ymax=None):
    f = sp.Abs(1-2*x)
    fn = sp.lambdify([x], f)
    psi, points= Lagrange_polynomials(x, N, [0, 1],
                                      point_distribution='Chebyshev')
    u = interpolation(f, psi, points)
    comparison_plot(f, u, Omega=[0, 1],
                    filename='Lagrange_interp_abs_Cheb_%d' % (N+1),
                    plot_title='Interpolation by Lagrange polynomials '\
                    'of degree %d' % N, ymin=ymin, ymax=ymax)
    print 'Interpolation points:', points

    # Make figures of Lagrange polynomials (psi)
    plt.figure()
    xcoor = np.linspace(0, 1, 1001)
    legends = []
    for i in (2, (N+1)/2+1):
        fn = sp.lambdify([x], psi[i])
        ycoor = fn(xcoor)
        plt.plot(xcoor, ycoor)
        legends.append(r'$\psi_%d$' % i)
        plt.hold('on')
    plt.legend(legends)
    plt.plot(points, [0]*len(points), 'ro')
    #if ymin is not None and ymax is not None:
    #    axis([xcoor[0], xcoor[-1], ymin, ymax])
    plt.savefig('Lagrange_basis_Cheb_%d.pdf' % (N+1))
    plt.savefig('Lagrange_basis_Cheb_%d.png' % (N+1))

def error_of_methods():
	"""
	Complie test with returning error
	"""
	import math
	import numpy as np
	import sympy
	from sympy import cos
	from sympy import sin
	from sympy.utilities.lambdify import lambdify
	
	w,x, L ,t= sympy.symbols("w x L t")
	pi = math.pi	

	u = lambda x,t : cos(w*t)*cos(pi*x/L)
	q = lambda x : 1 + ( x -L/2.0)**4

	def source_term(u,q):
		return sympy.simplify(u(x,t).diff(t, t) - (u(x,t).diff(x)*q(x)).diff(x) )

	w=1
	L = 2
	T=10
	Nx=40

	f =  sympy.lambdify((x,t), source_term(u,q)  ,'numpy')
	q = sympy.lambdify(x,q(x) ,'numpy')
	I = sympy.lambdify(x, u(x,t).subs(t,0),'numpy')
	u = sympy.lambdify((x,t), u(x,t),'numpy')

	print solve_wave_eqn_with_variable_velocity(q,f,I, L,T,Nx,u_exact=u,neumann=sum_approximation)[1]
	print solve_wave_eqn_with_variable_velocity(q,f,I, L,T,Nx,u_exact=u,neumann=centered_difference)[1]
	print solve_wave_eqn_with_variable_velocity(q,f,I, L,T,Nx,u_exact=u,neumann=centered_difference)[1]
	print solve_wave_eqn_with_variable_velocity(q,f,I, L,T,Nx,u_exact=u,neumann=shifted_domain)[1]

def comparison_plot(f, u, Omega, filename='tmp.pdf',
                    plot_title='', ymin=None, ymax=None,
                    u_legend='approximation'):
    """Compare f(x) and u(x) for x in Omega in a plot."""
    x = sm.Symbol('x')
    print 'f:', f

    f = sm.lambdify([x], f, modules="numpy")
    u = sm.lambdify([x], u, modules="numpy")
    if len(Omega) != 2:
        raise ValueError('Omega=%s must be an interval (2-list)' % str(Omega))
    # When doing symbolics, Omega can easily contain symbolic expressions,
    # assume .evalf() will work in that case to obtain numerical
    # expressions, which then must be converted to float before calling
    # linspace below
    if not isinstance(Omega[0], (int,float)):
        Omega[0] = float(Omega[0].evalf())
    if not isinstance(Omega[1], (int,float)):
        Omega[1] = float(Omega[1].evalf())

    resolution = 401  # no of points in plot
    xcoor = linspace(Omega[0], Omega[1], resolution)
    # Vectorized functions expressions does not work with
    # lambdify'ed functions without the modules="numpy"
    exact  = f(xcoor)
    approx = u(xcoor)
    plot(xcoor, approx, '-')
    hold('on')
    plot(xcoor, exact, '-')
    legend([u_legend, 'exact'])
    title(plot_title)
    xlabel('x')
    if ymin is not None and ymax is not None:
        axis([xcoor[0], xcoor[-1], ymin, ymax])
    savefig(filename)

 def fbenchmark(f, var=[Symbol('x')]):
     """
     Do some benchmarks with f using clambdify, lambdify and psyco.
     """
     global cf, pf, psyf
     start = time()
     cf = clambdify(var, f)
     print('compile time (including sympy overhead): %f s' % (
         time() - start))
     pf = lambdify(var, f, 'math')
     psyf = None
     psyco = import_module('psyco')
     if psyco:
         psyf = lambdify(var, f, 'math')
         psyco.bind(psyf)
     code = '''for x in (i/1000. for i in range(1000)):
     f(%s)''' % ('x,'*len(var)).rstrip(',')
     t1 = Timer(code, 'from __main__ import cf as f')
     t2 = Timer(code, 'from __main__ import pf as f')
     if psyf:
         t3 = Timer(code, 'from __main__ import psyf as f')
     else:
         t3 = None
     print('for x = (0, 1, 2, ..., 999)/1000')
     print('20 times in 3 runs')
     print('compiled:      %.4f %.4f %.4f' % tuple(t1.repeat(3, 20)))
     print('Python lambda: %.4f %.4f %.4f' % tuple(t2.repeat(3, 20)))
     if t3:
         print('Psyco lambda:  %.4f %.4f %.4f' % tuple(t3.repeat(3, 20)))

def derivative_matrix(deg, basis_functions, unroll=True):
    '''
    Matrix of \int_{-1}^{1}(L_i, L_j`) for polynomials spanned by deg+1 
    basis_functions.
    '''
    # Numerical quadrature, mas degree of integrand 2*deg - 1 -> deg
    xq, wq = np.polynomial.legendre.leggauss(deg)
    # Fast eval of basis and asis derivative
    x = Symbol('x')
    basis = [lambdify(x, f, 'numpy') for f in basis_functions(deg)]
    dbasis = [lambdify(x, f.diff(x, 1), 'numpy') for f in basis_functions(deg)]
    
    # Save calls to eval
    if unroll:
        V = np.zeros((len(basis), len(xq)))
        # Basis functions evaluated in quadrature points
        for row, f in enumerate(basis): V[row, :] = f(xq)

        dV = np.zeros((len(xq), len(basis)))
        # Basis derivatives evaluated in quadrature points
        for col, df in enumerate(dbasis): dV[:, col] = df(xq)
        
        # Integrate 
        C = (wq*V).dot(dV)

    # Here there are more calls to eval
    else:
        C = np.zeros((len(basis), len(basis)))
        for i, v in enumerate(basis):
            for j, du in enumerate(dbasis):
                C[i, j] = np.sum(wq*v(xq)*du(xq))
    
    return C

def comparison_plot(f, u, Omega, plotfile='tmp'):
    """Compare f(x,y) and u(x,y) for x,y in Omega in a plot."""
    x, y = sm.symbols('x y')

    f = sm.lambdify([x,y], f, modules="numpy")
    u = sm.lambdify([x,y], u, modules="numpy")
    # When doing symbolics, Omega can easily contain symbolic expressions,
    # assume .evalf() will work in that case to obtain numerical
    # expressions, which then must be converted to float before calling
    # linspace below
    for r in range(2):
        for s in range(2):
            if not isinstance(Omega[r][s], (int,float)):
                Omega[r][s] = float(Omega[r][s].evalf())

    resolution = 41  # no of points in plot
    xcoor = linspace(Omega[0][0], Omega[0][1], resolution)
    ycoor = linspace(Omega[1][0], Omega[1][1], resolution)
    xv, yv = ndgrid(xcoor, ycoor)
    # Vectorized functions expressions does not work with
    # lambdify'ed functions without the modules="numpy"
    exact  = f(xv, yv)
    approx = u(xv, yv)
    figure()
    surfc(xv, yv, exact, title='f(x,y)',
          colorbar=True, colormap=hot(), shading='flat')
    if plotfile:
        savefig('%s_f.pdf' % plotfile, color=True)
        savefig('%s_f.png' % plotfile)
    figure()
    surfc(xv, yv, approx, title='f(x,y)',
          colorbar=True, colormap=hot(), shading='flat')
    if plotfile:
        savefig('%s_u.pdf' % plotfile, color=True)
        savefig('%s_u.png' % plotfile)

def test_curly_matrix_symbol():
    # Issue #15009
    curlyv = sympy.MatrixSymbol("{v}", 2, 1)
    lam = lambdify(curlyv, curlyv)
    assert lam(1)==1
    lam = lambdify(curlyv, curlyv, dummify=True)
    assert lam(1)==1

def test_quadratic():
    """Verify a quadratic solution."""
    I = 1.2; V = 3; m = 2; b = 0.9
    s = lambda u: 4*u
    t = sym.Symbol('t')
    dt = 0.2
    T = 2

    q = 2  # arbitrary constant
    u_exact = I + V*t + q*t**2
    F = sym.lambdify(t, lhs_eq(t, m, b, s, u_exact, 'linear'))
    u_exact = sym.lambdify(t, u_exact, modules='numpy')
    #u1, t1 = solver(I, V, m, b, s, F, dt, T, 'linear')
    u1, t1 = solver_bwdamping(I, V, m, b, s, F, dt, T, 'linear')
    diff = np.abs(u_exact(t1) - u1).max()
    print diff
    tol = 1E-13
    #assert diff < tol

    # In the quadratic damping case, u_exact must be linear
    # in order to exactly recover this solution
    u_exact = I + V*t
    F = sym.lambdify(t, lhs_eq(t, m, b, s, u_exact, 'quadratic'))
    u_exact = sym.lambdify(t, u_exact, modules='numpy')
    #u2, t2 = solver(I, V, m, b, s, F, dt, T, 'quadratic')
    u2, t2 = solver_bwdamping(I, V, m, b, s, F, dt, T, 'quadratic')
    diff = np.abs(u_exact(t2) - u2).max()
    print diff
    assert diff < tol

    def make_numpy_fns_of_d1d2xw(self, dg_first, dg_second):

        args_x = self.xvar_tp1_sym + self.xvar_t_sym + self.xvar_tm1_sym

        args_w = self.wvar_tp1_sym + self.wvar_t_sym 

        args = args_x + args_w + self.param_sym_dict.values() 

        # args_values_x =   [x.subs(self.normal_xw_s_ss_values_d)
        #                  for x in args_x]               

        # args_values_w =   [x.subs(self.normal_xw_s_ss_values_d)
        #                  for x in args_w]               
        
        # args_values_p = [p.subs(self.par_to_values_dict)
        #                  for p in self.param_sym_dict.values()]

        # args_values = args_values_x + args_values_w + args_values_p

        dg_first_lam = sympy.lambdify(args, dg_first)
        dg_second_lam = sympy.lambdify(args, dg_second)

        self.fun_d_first_numpy = dg_first_lam
        self.fun_d_second_numpy = dg_second_lam

        return dg_first_lam, dg_second_lam

    def make_numpy_fns_of_d1d2xw(self, dg_first, dg_second):

        args = self.normal_x_s_d.values() + self.normal_w_s_d.values() + \
            self.param_sym_dict.values()

        args_values_x = [x.subs(self.normal_xw_s_ss_values_d)
                         for x in self.normal_x_s_d.values()]

        args_values_w = [w.subs(self.normal_xw_s_ss_values_d)
                         for w in self.normal_w_s_d.values()]

        args_values_p = [p.subs(self.par_to_values_dict)
                         for p in self.param_sym_dict.values()]

        args_values = args_values_x + args_values_w + args_values_p

#        print '\nargs for lambdify:', args
#        print '\nargs values for lambdify:', args_values
#        print '\nself.normal_xw_s_ss_values_d:', self.normal_xw_s_ss_values_d
#        print '\nself.par_to_values_dict:', self.par_to_values_dict

        dg_first_lam = sympy.lambdify(args, dg_first, dummify=False)
        dg_second_lam = sympy.lambdify(args, dg_second, dummify=False)

        return dg_first_lam, dg_second_lam, args_values

def make_integral(f_str):
    x          = sympy.Symbol('x')
    func_expr  = sympy.sympify(f_str)
    f          = sympy.lambdify(x, func_expr)
    int_expr   = sympy.integrate(func_expr)
    F          = sympy.lambdify(x, int_expr)
    return f, F

def make_derivative(f_str):
    x          = sympy.Symbol('x')
    func_expr  = sympy.sympify(f_str)
    f          = sympy.lambdify(x, func_expr)
    deriv_expr = sympy.diff(func_expr)
    df         = sympy.lambdify(x, deriv_expr)
    return f, df

def test_scipy_fns():
    if not scipy:
        skip("scipy not installed")

    single_arg_sympy_fns = [erf, erfc, factorial, gamma, loggamma, digamma]
    single_arg_scipy_fns = [scipy.special.erf, scipy.special.erfc,
        scipy.special.factorial, scipy.special.gamma, scipy.special.gammaln,
        scipy.special.psi]
    numpy.random.seed(0)
    for (sympy_fn, scipy_fn) in zip(single_arg_sympy_fns, single_arg_scipy_fns):
        test_values = 20 * numpy.random.rand(20)
        f = lambdify(x, sympy_fn(x), modules = "scipy")
        assert numpy.all(abs(f(test_values) - scipy_fn(test_values)) < 1e-15)

    double_arg_sympy_fns = [RisingFactorial, besselj, bessely, besseli,
        besselk]
    double_arg_scipy_fns = [scipy.special.poch, scipy.special.jn,
        scipy.special.yn, scipy.special.iv, scipy.special.kn]

    #suppress scipy warnings
    import warnings
    warnings.filterwarnings('ignore', '.*floating point number truncated*')

    for (sympy_fn, scipy_fn) in zip(double_arg_sympy_fns, double_arg_scipy_fns):
        for i in range(20):
            test_values = 20 * numpy.random.rand(2)
            f = lambdify((x,y), sympy_fn(x,y), modules = "scipy")
            assert abs(f(*test_values) - scipy_fn(*test_values)) < 1e-15

    def initialize(self, args):
        self.E   = args[0]
        self.A   = args[1]
        self.I   = args[2]
        self.r   = args[3]
        self.rho = args[4]
        self.l   = args[5]
        self.g   = args[6]

        q = sym.Matrix(sym.symarray('q',6))
        E, A, I, r, rho, l, g = sym.symbols('E A I r rho l g')
        theta = sym.Matrix(['theta_1','theta_2'])
        omega = sym.Matrix(['omega_1','omega_2'])
        
        # Load symbolic needed matricies and vectors
        M_sym      = pickle.load( open( "gebf-mass-matrix.dump",  "rb" ) )
        beta_sym   = pickle.load( open( "gebf-force-vector.dump", "rb" ) )
        Gamma1_sym = pickle.load( open( "gebf-1c-matrix.dump",    "rb" ) )
        Gamma2_sym = pickle.load( open( "gebf-2c-matrix.dump",    "rb" ) )

        # Create numeric function of needed matrix and vector quantities
        # this is MUCH MUCH faster than .subs()
        # .subs() is unusably slow !!
        self.M      = lambdify((E, A, I, r, rho, l, g, q, theta, omega),      M_sym, "numpy")
        self.beta   = lambdify((E, A, I, r, rho, l, g, q, theta, omega),   beta_sym, "numpy")
        self.Gamma1 = lambdify((E, A, I, r, rho, l, g, q, theta, omega), Gamma1_sym, "numpy")
        self.Gamma2 = lambdify((E, A, I, r, rho, l, g, q, theta, omega), Gamma2_sym, "numpy")