def __init__(self, nddistr=None, d=None):
        super(TwoVarsModel, self).__init__(nddistr, [d])
        self.eliminate_other([d])
        self.d = d
        self.vars = []
        self.symvars = []
        for var in nddistr.Vars: #self.free_rvs:
            self.vars.append(var)
            self.symvars.append(var.getSymname())
        #print "=====", self.vars
        #print self.symvars
        #print self.dep_rvs
        #print self.rv_to_equation
        self.symop = self.rv_to_equation[d]

        if len(self.vars) != 2:
            raise Exception("use it with two variables")
        x = self.symvars[0]
        y = self.symvars[1]
        z = sympy.Symbol("z")
        self.fun_alongx = eq_solve(self.symop, z, y)[0]
        self.fun_alongy = eq_solve(self.symop, z, x)[0]

        self.lfun_alongx = my_lambdify([x, z], self.fun_alongx, "numpy")
        self.lfun_alongy = my_lambdify([y, z], self.fun_alongy, "numpy")
        self.Jx = 1 * sympy.diff(self.fun_alongx, z)
        #print "Jx=", self.Jx
        #print "fun_alongx=", self.fun_alongx
        self.Jy = 1 * sympy.diff(self.fun_alongy, z)
        self.lJx = my_lambdify([x, z], self.Jx, "numpy")
        self.lJy = my_lambdify([y, z], self.Jy, "numpy")
        self.z = z

def test_atan2():
    assert atan2(0, 0) == S.NaN
    assert atan2(0, 1) == 0
    assert atan2(1, 1) == pi/4
    assert atan2(1, 0) == pi/2
    assert atan2(1, -1) == 3*pi/4
    assert atan2(0, -1) == pi
    assert atan2(-1, -1) == -3*pi/4
    assert atan2(-1, 0) == -pi/2
    assert atan2(-1, 1) == -pi/4

    u = Symbol("u", positive=True)
    assert atan2(0, u) == 0
    u = Symbol("u", negative=True)
    assert atan2(0, u) == pi

    assert atan2(y, oo) ==  0
    assert atan2(y, -oo)==  2*pi*Heaviside(re(y)) - pi

    assert atan2(y, x).rewrite(log) == -I*log((x + I*y)/sqrt(x**2 + y**2))
    assert atan2(y, x).rewrite(atan) == 2*atan(y/(x + sqrt(x**2 + y**2)))

    assert diff(atan2(y, x), x) == -y/(x**2 + y**2)
    assert diff(atan2(y, x), y) == x/(x**2 + y**2)

    assert simplify(diff(atan2(y, x).rewrite(log), x)) == -y/(x**2 + y**2)
    assert simplify(diff(atan2(y, x).rewrite(log), y)) ==  x/(x**2 + y**2)

    assert isinstance(atan2(2, 3*I).n(), atan2)

def Simple_manifold_with_scalar_function_derivative():
    Print_Function()
    coords = (x,y,z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3',metric='[1,1,1]',coords=coords)
    # Define surface
    mfvar = (u,v) = symbols('u v')
    X = u*e1+v*e2+(u**2+v**2)*e3
    print '\\f{X}{u,v} =',X
    MF = Manifold(X,mfvar)
    (eu,ev) = MF.Basis()
    # Define field on the surface.
    g = (v+1)*log(u)

    print '\\f{g}{u,v} =',g

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g,u) + (MF.rbasis[1]/MF.E_sq)*diff(g,v)
    print '\\eval{\\nabla g}{u=1,v=0} =', VectorDerivative.subs({u:1,v:0})

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print '\\eval{\\f{Grad}{g}}{u=1,v=0} =', dg.subs({u:1,v:0})
    dg = MF.grad*g
    print '\\eval{\\nabla g}{u=1,v=0} =', dg.subs({u:1,v:0})
    return

    def __init__(self):
        self.num_lines, self.num_marks = 5, 2
        self._qp = {}

        function_type = random.choice(['sqrt', 'quadratic', 'product'])

        if function_type == 'sqrt':
            outer_function = all_functions.request_linear(difficulty=2).equation
            inner_function = all_functions.request_linear(difficulty=1).equation
            inner_function = inner_function.replace(lambda expr: expr.is_Symbol, lambda expr: sympy.sqrt(expr))

            self._qp['equation'] = outer_function.replace(x, inner_function)
            self._qp['derivative'] = sympy.diff(self._qp['equation'])

        elif function_type == 'quadratic':
            power_two_coeff = not_named_yet.randint_no_zero(-3, 3)
            power_one_coeff = not_named_yet.randint_no_zero(-5, 5)
            inner_function = power_two_coeff * x ** 2 + power_one_coeff * x
            index = random.randint(3, 5)

            self._qp['equation'] = inner_function ** index
            self._qp['derivative'] = sympy.diff(self._qp['equation'])

        elif function_type == 'product':
            left_function = x ** random.randint(1, 3)
            right_outer_function = random.choice([sympy.sin, sympy.cos, sympy.log, sympy.exp])
            right_inner_function = not_named_yet.randint_no_zero(-3, 3) * x

            self._qp['equation'] = left_function * right_outer_function(right_inner_function)
            self._qp['derivative'] = sympy.diff(self._qp['equation'])

def test_derivative_evaluate():
    assert Derivative(sin(x), x) != diff(sin(x), x)
    assert Derivative(sin(x), x).doit() == diff(sin(x), x)

    assert Derivative(Derivative(f(x), x), x) == diff(f(x), x, x)
    assert Derivative(sin(x), x, 0) == sin(x)
    assert Derivative(sin(x), (x, y), (x, -y)) == sin(x)

 def grad(self, func):
     """
     Calculate the gradient of 'func'.
     """
     return Matrix([diff(func, self.xs[0]) / self.h[0], \
                    diff(func, self.xs[1]) / self.h[1], \
                    diff(func, self.xs[2]) / self.h[2]])

def test_gegenbauer():
    n = Symbol("n")
    a = Symbol("a")

    assert gegenbauer(0, a, x) == 1
    assert gegenbauer(1, a, x) == 2*a*x
    assert gegenbauer(2, a, x) == -a + x**2*(2*a**2 + 2*a)
    assert gegenbauer(3, a, x) == \
        x**3*(4*a**3/3 + 4*a**2 + 8*a/3) + x*(-2*a**2 - 2*a)

    assert gegenbauer(-1, a, x) == 0
    assert gegenbauer(n, S(1)/2, x) == legendre(n, x)
    assert gegenbauer(n, 1, x) == chebyshevu(n, x)
    assert gegenbauer(n, -1, x) == 0

    X = gegenbauer(n, a, x)
    assert isinstance(X, gegenbauer)

    assert gegenbauer(n, a, -x) == (-1)**n*gegenbauer(n, a, x)
    assert gegenbauer(n, a, 0) == 2**n*sqrt(pi) * \
        gamma(a + n/2)/(gamma(a)*gamma(-n/2 + S(1)/2)*gamma(n + 1))
    assert gegenbauer(n, a, 1) == gamma(2*a + n)/(gamma(2*a)*gamma(n + 1))

    assert gegenbauer(n, Rational(3, 4), -1) == zoo

    m = Symbol("m", positive=True)
    assert gegenbauer(m, a, oo) == oo*RisingFactorial(a, m)

    assert conjugate(gegenbauer(n, a, x)) == gegenbauer(n, conjugate(a), conjugate(x))

    assert diff(gegenbauer(n, a, x), n) == Derivative(gegenbauer(n, a, x), n)
    assert diff(gegenbauer(n, a, x), x) == 2*a*gegenbauer(n - 1, a + 1, x)

def Delp2(f):
    """ Laplacian in X-Z
    """
    d2fdx2 = diff(f, x, 2)
    d2fdz2 = diff(f, z, 2)

    return d2fdx2 + d2fdz2

def ItoFormula_2(function,a,b,dt):
	func = compile(function,"",'eval')
	a=float(a)
	b=float(b)
	n=int((b-a)/dt)	
	WP = 0
	t = a
	for i in range(int(a/dt)):
		WP = WP+np.random.normal(0,math.sqrt(dt))
	part_1 = eval(func)
	f = sym.sympify(function)
	WP = sym.Symbol('WP')
	df = str(sym.diff(f,WP))
	ddf = str(sym.diff(df,WP))
	part_2 = 0
	part_3 = 0
	result_2 = 0
	WP = 0
	for i in range(0,n):
		t=dt*((a/dt)+i)
		WP = WP+np.random.normal(0,math.sqrt(dt))
		k2 = eval(df,valueDic)
		k3 = eval(ddf,valueDic)
		part_2 = part_2 + k2*(WT[i+1]-WT[i])
		part_3 = part_3 + k3*dt
	result_2 = part_1 + part_2 + part_3
	print 'ito formula result: '+str(result_2)

def eleq(Lagrangian, Friction = 0, t = Symbol('t')):
    """
    Returns Euler-Lagrange equations of the lagrangian system.

    Examples
    ========

    >>> from sympy import *
    >>> t, k = symbols('t k')
    >>> x = symbols('x', cls=Function)
    >>> eleq(diff(x(t),t)**2/2 - k*x(t)**2/2)
    {a_x: -k*x}

    >>> a = symbols('a')
    >>> eleq(diff(x(t),t)**2/2 - k*x(t)**2/2, a*diff(x(t),t)**2/2)
    {a_x: -*a*v_x - k*x}
    """
    Lagrangian = simplify(Lagrangian)
    var_list = [list(x.atoms(Function))[0] for x in Lagrangian.atoms(Derivative)]
    nvar = len(var_list)
    ecu_list = [ diff(Lagrangian, variable) - diff(Lagrangian, diff(variable,t), t) -  diff(Friction, diff(variable,t)) for variable in var_list ]
    str_list = [ str(variable).replace("("+str(t)+")","") for variable in var_list ]
    a_subs = {diff(var_list[i],t,2): Symbol('a_' + str_list[i]) for i in range(nvar)}
    v_subs = {diff(var_list[i],t): Symbol('v_' + str_list[i]) for i in range(nvar)}
    x_subs = {var_list[i]: Symbol(str_list[i]) for i in range(nvar)}
    for i in range(nvar):
        if hasattr(ecu_list[i], "subs"):
            ecu_list[i] = ecu_list[i].subs(a_subs).subs(v_subs).subs(x_subs)
    a_list = sorted(list(a_subs.values()), key = str)
    return solveswc(ecu_list, a_list)

    def _derive_dw_dpsi():
        """Derive the slope dw_dpsi with sympy
        
        Returns
        -------
        dw_dpsi1 : sympy expression
            Derivative including correction factor.
        dw_dpsi2 : sympy expression
            Derivative with correction factor equal to one.
         """
        psi, a, n, m, psir, ws = sympy.symbols('psi, a, n, m, psir, ws')
        psi, ws, a, b, wr, s1, psir = sympy.symbols('psi, ws, a, b, wr, s1, psir')
        from sympy import log, exp        
        C1 = (1 - log(1 + psi / psir) / 
                          log(1.0 + 1.0e6 / psir))

        C2 = 1      
        
        l10 = sympy.symbols('l10')
        
        w = C1 * ((ws - s1 * log(psi)/l10 - wr) * a / (psi**b + a) + wr)
        dw_dpsi1 = sympy.diff(w, psi)                    
        
        w = C2 * ((ws - s1 * log(psi)/l10 - wr) * a / (psi**b + a) + wr)
        dw_dpsi2 = sympy.diff(w, psi)
        return  dw_dpsi1, dw_dpsi2

def hamiltonian(Lagrangian, t = Symbol('t'), delta = False):
    """
    Returns the Hamiltonian of the Lagrangian.

    Examples
    ========

    >>> from sympy import *
    >>> t, k = symbols('t k')
    >>> x = symbols('x', cls=Function)
    >>> hamiltonian(diff(x(t),t)**2/2 - k*x(t)**2/2)
    k*x**2/2 + v_x**2/2
    """
    Lagrangian = simplify(Lagrangian)
    var_list = [list(x.atoms(Function))[0] for x in Lagrangian.atoms(Derivative)]
    nvar = len(var_list)
    # New variables.
    str_list = [ str(variable).replace("("+str(t)+")","") for variable in var_list ]
    v_subs = {diff(var_list[i],t): Symbol('v_' + str_list[i]) for i in range(nvar)}
    x_subs = {var_list[i]: Symbol(str_list[i]) for i in range(nvar)}
    # Hamiltonian calculus.
    dxdLv = 0
    for variable in var_list:
        dxdLv += diff(variable,t)*diff(Lagrangian, diff(variable,t))
    result = simplify((dxdLv - Lagrangian).subs(v_subs).subs(x_subs))
    if delta:
        v0_subs = {Symbol('v_' + str_list[i]): Symbol('v_' + str_list[i] + "0") for i in range(nvar)}
        x0_subs = {Symbol(str_list[i]): Symbol(str_list[i] + "0") for i in range(nvar)}
        return result - result.subs(v0_subs).subs(x0_subs)
    else:
        return result

def test_legendre():
    raises(ValueError, lambda: legendre(-1, x))
    assert legendre(0, x) == 1
    assert legendre(1, x) == x
    assert legendre(2, x) == ((3*x**2 - 1)/2).expand()
    assert legendre(3, x) == ((5*x**3 - 3*x)/2).expand()
    assert legendre(4, x) == ((35*x**4 - 30*x**2 + 3)/8).expand()
    assert legendre(5, x) == ((63*x**5 - 70*x**3 + 15*x)/8).expand()
    assert legendre(6, x) == ((231*x**6 - 315*x**4 + 105*x**2 - 5)/16).expand()

    assert legendre(10, -1) == 1
    assert legendre(11, -1) == -1
    assert legendre(10, 1) == 1
    assert legendre(11, 1) == 1
    assert legendre(10, 0) != 0
    assert legendre(11, 0) == 0

    assert roots(legendre(4, x), x) == {
        sqrt(Rational(3, 7) - Rational(2, 35)*sqrt(30)): 1,
        -sqrt(Rational(3, 7) - Rational(2, 35)*sqrt(30)): 1,
        sqrt(Rational(3, 7) + Rational(2, 35)*sqrt(30)): 1,
        -sqrt(Rational(3, 7) + Rational(2, 35)*sqrt(30)): 1,
    }

    n = Symbol("n")

    X = legendre(n, x)
    assert isinstance(X, legendre)

    assert legendre(-n, x) == legendre(n - 1, x)
    assert legendre(n, -x) == (-1)**n*legendre(n, x)
    assert diff(legendre(n, x), x) == \
        n*(x*legendre(n, x) - legendre(n - 1, x))/(x**2 - 1)
    assert diff(legendre(n, x), n) == Derivative(legendre(n, x), n)

def make_qdiff_to_q012(x_names):
    qdiff_to_q012 = {}
    qdiff_0 = {}
    qdiff_1 = {}
    qdiff_2 = {}

    for date_str in ['tp1', 't', 'tm1']:
        dated_names = [name + date_str for name in x_names]
        dated_q_fun_sym = [sympy.Function(x)(q) for x in dated_names]

        dated_names_0 = [name + '_0' for name in dated_names]
        dated_names_0_sym = [sympy.Symbol(x) for x in dated_names_0]
        qdiff_0.update(dict(zip(dated_q_fun_sym, dated_names_0_sym)))

        dated_qdiffs_1 = [sympy.diff(x, q, 1) for x in dated_q_fun_sym]
        dated_names_1 = [name + '_1' for name in dated_names]
        dated_names_1_sym = [sympy.Symbol(x) for x in dated_names_1]
        qdiff_1.update(dict(zip(dated_qdiffs_1, dated_names_1_sym)))

        dated_qdiffs_2 = [sympy.diff(x, q, 2) for x in dated_q_fun_sym]
        dated_names_2 = [name + '_2' for name in dated_names]
        dated_names_2_sym = [sympy.Symbol(x) for x in dated_names_2]
        qdiff_2.update(dict(zip(dated_qdiffs_2, dated_names_2_sym)))

    qdiff_to_q012.update(qdiff_0)
    qdiff_to_q012.update(qdiff_1)
    qdiff_to_q012.update(qdiff_2)

    return qdiff_to_q012

    def __call__(self, scalar_field):
        """
        Represents the gradient of the given scalar field.

        Parameters
        ==========

        scalar_field : SymPy expression
            The scalar field to calculate the gradient of.

        Examples
        ========

        >>> from sympy.vector import CoordSysCartesian
        >>> C = CoordSysCartesian('C')
        >>> C.delop(C.x*C.y*C.z)
        C.y*C.z*C.i + C.x*C.z*C.j + C.x*C.y*C.k

        """

        scalar_field = express(scalar_field, self.system,
                               variables = True)
        vx = diff(scalar_field, self._x)
        vy = diff(scalar_field, self._y)
        vz = diff(scalar_field, self._z)

        return vx*self._i + vy*self._j + vz*self._k

def differentiate(C,a, b):
    diff1 = [0]*N
    diff1_anal = [0]*N
    diff2 = [0]*N
    diff2_anal = [0]*N
    i = arange(N)
    # Chebyshev nodes on [-1,1]
    x = -cos(pi*(2.0*i + 1.0)/(2.0*(N-1) + 2.0)) 

    #chebdiff1 = pow(2,1-i)*i*sin(i*arccos(x))/pow(1-x*x,0.5)
    #chebdiff2 = pow(2,1-i)*i*( x*sin(i*arccos(x))/pow(1-x*x,1.5) - i*cos(i*arccos(x))/(1-x*x))

    # scale to physical interval [a,b]
    xp = affineTransform(x, a, b)
    
    c = 2.0/(b-a)
    for i in range(0,N):
        for j in range(0, N):
            tmp = pow(2,1-j)*j*sin(j*arccos(x[i]))/pow(1-x[i]*x[i],0.5)
            diff1[i] += c*C[j]*tmp
            tmp = pow(2,1-j)*j*( x[i]*sin(j*arccos(x[i]))/pow(1-x[i]*x[i],1.5) - j*cos(j*arccos(x[i]))/(1-x[i]*x[i]))
            diff2[i] += c*c*C[j]*tmp
        diff1_anal[i] = (sympy.diff(sympy.exp(-sympy.abc.x*sympy.abc.x*beta), sympy.abc.x, 1)).subs(sympy.abc.x,xp[i])
        diff2_anal[i] = (sympy.diff(sympy.exp(-sympy.abc.x*sympy.abc.x*beta), sympy.abc.x, 2)).subs(sympy.abc.x,xp[i])
        
    # plot the Chebyshev nodes and the interpolation
    plot(xp, diff1, 'r--', label = "chebyshev derivative ord 1")
    plot(xp, diff1_anal, 'g.', label = "analytic derivative ord 1")
    plot(xp, diff2, 'y--', label = "chebyshev derivative ord 2")
    plot(xp, diff2_anal, 'b*', label = "analytic derivative ord 2")
    legend(loc="best")
    grid(True)
    show()

 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 update_integral_data(u_str, v_str, cx_str, cy_str):
    #string parsing
    u_fun, _ = my_bokeh_utils.string_to_function_parser(u_str, ['x', 'y'])
    v_fun, _ = my_bokeh_utils.string_to_function_parser(v_str, ['x', 'y'])
    cx_fun, cx_sym = my_bokeh_utils.string_to_function_parser(cx_str, ['t'])
    cy_fun, cy_sym = my_bokeh_utils.string_to_function_parser(cy_str, ['t'])

    from sympy import diff

    dcx_sym = diff(cx_sym,'t')
    dcy_sym = diff(cy_sym,'t')

    dcx_fun = my_bokeh_utils.sym_to_function_parser(dcx_sym, 't')
    dcy_fun = my_bokeh_utils.sym_to_function_parser(dcy_sym, 't')

    t = np.linspace(curveintegral_settings.parameter_min,curveintegral_settings.parameter_max)

    f_I = lambda xx, tt: (u_fun(cx_fun(tt),cy_fun(tt)) * dcx_fun(tt) + v_fun(cx_fun(tt),cy_fun(tt)) * dcy_fun(tt))

    integrand = f_I(None,t)

    from scipy.integrate import odeint

    integral = odeint(f_I,0,t)

    source_integral.data = dict(t=t.tolist(),
                                integral=integral.tolist(),
                                integrand=integrand.tolist())

def case0(f, N=3):
    B = 1 - x ** 3
    dBdx = sm.diff(B, x)

    # Compute basis functions and their derivatives
    phi = {0: [x ** (i + 1) * (1 - x) for i in range(N + 1)]}
    phi[1] = [sm.diff(phi_i, x) for phi_i in phi[0]]

    def integrand_lhs(phi, i, j):
        return phi[1][i] * phi[1][j]

    def integrand_rhs(phi, i):
        return f * phi[0][i] - dBdx * phi[1][i]

    Omega = [0, 1]

    u_bar = solve(integrand_lhs, integrand_rhs, phi, Omega, verbose=True, numint=False)
    u = B + u_bar
    print "solution u:", sm.simplify(sm.expand(u))

    # Calculate analytical solution

    # Solve -u''=f by integrating f twice
    f1 = sm.integrate(f, x)
    f2 = sm.integrate(f1, x)
    # Add integration constants
    C1, C2 = sm.symbols("C1 C2")
    u_e = -f2 + C1 * x + C2
    # Find C1 and C2 from the boundary conditions u(0)=0, u(1)=1
    s = sm.solve([u_e.subs(x, 0) - 1, u_e.subs(x, 1) - 0], [C1, C2])
    # Form the exact solution
    u_e = -f2 + s[C1] * x + s[C2]
    print "analytical solution:", u_e
    # print 'error:', u - u_e  # many terms - which cancel
    print "error:", sm.expand(u - u_e)

def update_parameter_data(cx_str, cy_str, t):
    # string parsing
    cx_fun, cx_sym = my_bokeh_utils.string_to_function_parser(cx_str, ['t'])
    cy_fun, cy_sym = my_bokeh_utils.string_to_function_parser(cy_str, ['t'])

    from sympy import diff

    dcx_sym = diff(cx_sym, 't')
    dcy_sym = diff(cy_sym, 't')

    dcx_fun = my_bokeh_utils.sym_to_function_parser(dcx_sym, 't')
    dcy_fun = my_bokeh_utils.sym_to_function_parser(dcy_sym, 't')

    # crating sample
    x_val = cx_fun(t)
    y_val = cy_fun(t)

    dx_val = dcx_fun(t)
    dy_val = dcy_fun(t)
    xx, hx = np.linspace(source_view.data['x_start'][0], source_view.data['x_end'][0], curveintegral_settings.n_sample,
                         retstep=True)
    ssdict, spdict, _ = my_bokeh_utils.quiver_to_data(x=np.array(x_val), y=np.array(y_val), u=np.array(dx_val), v=np.array(dy_val), h=2*hx,
                                                           do_normalization=True, fix_at_middle=False)

    # save data
    source_param.data = dict(x=[x_val], y=[y_val], t=[t], x0=ssdict['x0'], y0=ssdict['y0'],
                             x1=ssdict['x1'], y1=ssdict['y1'], xs=spdict['xs'], ys=spdict['ys'])

    print "curve point was updated with t=%f" % (t)