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")

def velocity_field(psi): #takes a symbolic function and returns two lambda functions
#to evaluate the derivatives in both x and y.
   global w
   if velocity_components:
      u = lambdify((x,y), eval(x_velocity), modules='numpy')
      v = lambdify((x,y), eval(y_velocity), modules='numpy')
   else:
      if is_complex_potential:
         print "Complex potential, w(z) given"
         #define u, v symbolically as the imaginary part of the derivatives
         u = lambdify((x, y), sympy.im(psi.diff(y)), modules='numpy')
         v = lambdify((x, y), -sympy.im(psi.diff(x)), modules='numpy')
      else:
         #define u,v as the derivatives 
         print "Stream function, psi given"
         u = sympy.lambdify((x, y), psi.diff(y), 'numpy')
         v = sympy.lambdify((x, y), -psi.diff(x), 'numpy')
   if (branch_cuts): # If it's indicated that there are branch cuts in the mapping,
                      # then we need to return vectorized numpy functions to evaluate
                      # everything numerically, instead of symbolically 
                      # This of course results in a SIGNIFICANT time increase
                      #   (I don't know how to handle more than the primitive root
                      #   (symbolically in Sympy
      return np.vectorize(u),np.vectorize(v)
   else:
       # If there are no branch cuts, then return the symbolic lambda functions (MUCH faster)
      return u,v

    def calculate_path(self):

        f, g, t, dt, s = (self.f, self.g, self.t, self.delta, self.s)
        der_f, der_g = [diff(coord, t) for coord in f, g]
        tanf, tang = [lambdify((t, s), coord + der * s, numpy) for coord, der in ((f, der_f), (g, der_g))]

        coords = lambdify(t, (f, g), numpy)
        print der_f.subs(t, t + dt)
        get_control_param = lambdify(t, (g - g.subs(t, t + dt)) / (der_g.subs(t, t + dt) - der_g), numpy)

        def control_coords(time):
            control_param = get_control_param(time)
            return (tangent(time, control_param) for tangent in (tanf, tang))

        x_o, y_o = coords(self.min)
        path = "M%f, %f" % (x_o, y_o)

        def append_curve(path, point, control_point):
            return "%s Q%f, %f %f, %f" % (path, point.x, point.y, control_point.x, control_point.y)

        for time in numpy.arange(self.min + self.delta, self.max - self.delta, self.delta):
            x, y = coords(time)
            x_c, y_c = control_coords(time)
            path = append_curve(path, Point(x, y), Point(x_c, y_c))

        x_f, y_f = coords(self.max)
        path = "%s T%f, %f" % (path, x_f, y_f)

        self.path = path

    def __init__(self, s=1, h1=0.2, h2=0.6, h3=0.2, v1=1, v2=1, v3=1,
                 base=None):
        """Inputs: s  - the stratification parameter
                   h1 - depth of layer 1
                   h2 - depth of layer 2
                   h3 - depth of layer 3
                   v1 - speed relation for layer 1
                   v2 - speed relation for layer 2
                   v3 - speed relation for layer 3.
                   base - optional precalculated base equation set
        """
        # Create functions of (a, b) for these given parameters
        self.F = self.F(v1=v1, v2=v2, v3=v3, h1=h1, h2=h2, h3=h3, s=s)
        self.G = self.G(v1=v1, v2=v2, v3=v3, h1=h1, h2=h2, h3=h3)
        self.f = sp.lambdify((a, b), self.F)
        self.g = sp.lambdify((a, b), self.G)

        self.H = h1, h2, h3
        self.V = v1, v2, v3

        self.base = base or LambBase(s=s, h1=h1, h2=h2, h3=h3)

        # Sample grid for finding F(a, b) = 0
        # resolution of rough zero search (np.arange)
        self.resolution = 0.01

        # bounds on physical solutions
        a_min, a_max = (-h1, 1 - h1)
        b_min, b_max = (-(h1 + h2), h3)

        s = 2 * self.resolution  # extra edge to consider
        x = np.arange(a_min - s, a_max + s, self.resolution)
        y = np.arange(b_min - s, b_max + s, self.resolution)

        self.AB = np.meshgrid(x, y)

 def __init__(self,
              fi=lambda t, theta: log(t),
              fi_inv=None, #lambda t, theta:(-sympy.log(t)) ** theta,
              theta=2,
              marginals=None):
     self.theta = float(theta)#Symbol('theta')
     self.t = Symbol('t')
     self.s = Symbol('s')
     self.d = len(marginals)
     self.fi_ = fi
     self.fi_inv_ = fi_inv
     self.sym_fi = fi(self.t, self.theta)
     self.sym_fi_deriv = sympy.diff(self.sym_fi, self.t)
     if fi_inv is  None:
         self.sym_fi_inv = sympy.solve(self.sym_fi - self.s, self.t)[0]
     else:
         self.sym_fi_inv = fi_inv(self.s, self.theta)
     self.sym_fi_inv_nth_deriv = sympy.diff(self.sym_fi_inv, self.s, self.d)
     #self.debug_info()
     super(ArchimedeanSymbolicCopula, self).__init__(fi=sympy.lambdify(self.t, self.sym_fi, "numpy"),
                                                     fi_deriv=sympy.lambdify(self.t, self.sym_fi_deriv, "numpy"),
                                                     fi_inv=sympy.lambdify(self.s, self.sym_fi_inv, "numpy"),
                                                     fi_inv_nth_deriv=sympy.lambdify(self.s, self.sym_fi_inv_nth_deriv, "numpy"),
                                                     marginals=marginals)
     vars = self.symVars
     si = 0
     for i in range(self.d):
         si += self.fi_(vars[i], self.theta)
     self.sym_C = self.fi_inv_(si, self.theta)

    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
        try:
            import psyco

            psyf = lambdify(var, f, "math")
            psyco.bind(psyf)
        except ImportError:
            pass
        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))