Python sympy.N Examples

def _scheme_from_rc_mpmath(alpha, beta):
    # Create vector cut of the first value of beta
    n = len(alpha)
    b = mp.zeros(n, 1)
    for i in range(n - 1):
        b[i] = mp.sqrt(beta[i + 1])

    z = mp.zeros(1, n)
    z[0, 0] = 1
    d = mp.matrix(alpha)
    tridiag_eigen(mp, d, b, z)

    # nx1 matrix -> list of mpf
    x = numpy.array([mp.mpf(sympy.N(xx, mp.dps)) for xx in d])
    w = numpy.array([mp.mpf(sympy.N(beta[0], mp.dps)) * mp.power(ww, 2) for ww in z])
    return x, w 

def event_bounds_expressions(self, event_bounds_exp):
        if hasattr(self, 'output_equations'):
            assert len(event_bounds_exp)+1 == self.output_equations.shape[0]
        if hasattr(self, 'output_equations_functions'):
            assert len(event_bounds_exp)+1 == \
                self.output_equations_functions.size
        if getattr(self, 'state_equations', None) is not None:
            assert len(event_bounds_exp)+1 == self.state_equations.shape[0]
        if getattr(self, 'state_equations_functions', None) is not None:
            assert len(event_bounds_exp)+1 == \
                self.state_equations_functions.size
        self._event_bounds_expressions = event_bounds_exp
        self.event_bounds = np.array(
            [sp.N(bound, subs=self.constants_values)
             for bound in event_bounds_exp],
            dtype=np.float_
        ) 

def Nga(x, prec=5):
    """
    Like :func:`sympy.N`, but also works on multivectors

    For multivectors with coefficients that contain floating point numbers, this
    rounds all these numbers to a precision of ``prec`` and returns the rounded
    multivector.
    """
    if isinstance(x, Mv):
        return Mv(Nsympy(x.obj, prec), ga=x.Ga)
    else:
        return Nsympy(x, prec) 

def test_Wigner3j_values():
    from sympy import N
    from sympy.physics.wigner import wigner_3j
    from spherical_functions import Wigner3j
    j_max = 8
    for j1 in range(j_max+1):
        for j2 in range(j1, j_max+1):
            for j3 in range(j2, j_max+1):
                for m1 in range(-j1, j1 + 1):
                    for m2 in range(-j2, j2 + 1):
                        m3 = -m1-m2
                        if j3 >= abs(m3):
                            sf_3j = Wigner3j(j1, j2, j3, m1, m2, m3)
                            sy_3j = N(wigner_3j(j1, j2, j3, m1, m2, m3))
                            assert abs(sf_3j - sy_3j) < precision_Wigner3j 

def wigner_d_naive(l, m, n, beta):
    """
    Numerically naive implementation of the Wigner-d function.
    This is useful for checking the correctness of other implementations.

    :param l: the degree of the Wigner-d function. l >= 0
    :param m: the order of the Wigner-d function. -l <= m <= l
    :param n: the order of the Wigner-d function. -l <= n <= l
    :param beta: the argument. 0 <= beta <= pi
    :return: d^l_mn(beta) in the TODO: what basis? complex, quantum(?), centered, cs(?)
    """
    from scipy.special import eval_jacobi
    try:
        from scipy.misc import factorial
    except:
        from scipy.special import factorial

    from sympy.functions.special.polynomials import jacobi, jacobi_normalized
    from sympy.abc import j, a, b, x
    from sympy import N
    #jfun = jacobi_normalized(j, a, b, x)
    jfun = jacobi(j, a, b, x)
    # eval_jacobi = lambda q, r, p, o: float(jfun.eval(int(q), int(r), int(p), float(o)))
    # eval_jacobi = lambda q, r, p, o: float(N(jfun, int(q), int(r), int(p), float(o)))
    eval_jacobi = lambda q, r, p, o: float(jfun.subs({j:int(q), a:int(r), b:int(p), x:float(o)}))

    mu = np.abs(m - n)
    nu = np.abs(m + n)
    s = l - (mu + nu) / 2
    xi = 1 if n >= m else (-1) ** (n - m)

    # print(s, mu, nu, np.cos(beta), type(s), type(mu), type(nu), type(np.cos(beta)))
    jac = eval_jacobi(s, mu, nu, np.cos(beta))
    z = np.sqrt((factorial(s) * factorial(s + mu + nu)) / (factorial(s + mu) * factorial(s + nu)))

    # print(l, m, n, beta, np.isfinite(mu), np.isfinite(nu), np.isfinite(s), np.isfinite(xi), np.isfinite(jac), np.isfinite(z))
    assert np.isfinite(mu) and np.isfinite(nu) and np.isfinite(s) and np.isfinite(xi) and np.isfinite(jac) and np.isfinite(z)
    assert np.isfinite(xi * z * np.sin(beta / 2) ** mu * np.cos(beta / 2) ** nu * jac)
    return xi * z * np.sin(beta / 2) ** mu * np.cos(beta / 2) ** nu * jac 

def test_xk(k):
    n = 10

    moments = quadpy.tools.integrate(
        lambda x: [x ** (i + k) for i in range(2 * n)], -1, +1
    )

    alpha, beta = quadpy.tools.chebyshev(moments)

    assert (alpha == 0).all()
    assert beta[0] == moments[0]
    assert beta[1] == sympy.S(k + 1) / (k + 3)
    assert beta[2] == sympy.S(4) / ((k + 5) * (k + 3))
    quadpy.tools.scheme_from_rc(
        numpy.array([sympy.N(a) for a in alpha], dtype=float),
        numpy.array([sympy.N(b) for b in beta], dtype=float),
        mode="numpy",
    )

    def leg_polys(x):
        return orthopy.line_segment.tree_legendre(x, 19, "monic", symbolic=True)

    moments = quadpy.tools.integrate(
        lambda x: [x ** k * leg_poly for leg_poly in leg_polys(x)], -1, +1
    )

    _, _, a, b = orthopy.line_segment.recurrence_coefficients.legendre(
        2 * n, "monic", symbolic=True
    )

    alpha, beta = quadpy.tools.chebyshev_modified(moments, a, b)

    assert (alpha == 0).all()
    assert beta[0] == moments[0]
    assert beta[1] == sympy.S(k + 1) / (k + 3)
    assert beta[2] == sympy.S(4) / ((k + 5) * (k + 3))
    points, weights = quadpy.tools.scheme_from_rc(
        numpy.array([sympy.N(a) for a in alpha], dtype=float),
        numpy.array([sympy.N(b) for b in beta], dtype=float),
        mode="numpy",
    ) 

def get_generating_function(f, T=None):
    """
    Get the generating function

    :param f: function addr
    :param T: edges matrix of the function f
    :return: the generating function
    """

    if not T:
        edges = [(a[0].addr, a[1].addr) for a in f.graph.edges()]
        N = sorted([x for x in f.block_addrs])
        T = []
        for b1 in N:
            T.append([])
            for b2 in N:
                if b1 == b2 or (b1, b2) in edges:
                    T[-1].append(1)
                else:
                    T[-1].append(0)
    else:
        N = T[0]

    T = sympy.Matrix(T)
    z = sympy.var('z')
    I = sympy.eye(len(N))

    tmp = I - z * T
    tmp.row_del(len(N) - 1)
    tmp.col_del(0)
    det_num = tmp.det()
    det_den = (I - z * T).det()
    quot = det_num / det_den
    g_z = ((-1) ** (len(N) + 1)) * quot
    return g_z 

def get_path_n(f, T=None):
    # TODO cire paper here
    """
    Get the formula to estimate the number of path of a function, given its longest path.

    :param f: function addr
    :param T: edges matrix of the function f
    :return: formula to estimate the number of paths
    """

    g_z = get_generating_function(f, T)
    expr = g_z.as_numer_denom()[1]
    rs = sympy.roots(expr)
    D = len(set(rs.keys()))  # number of distinct roots
    d = sum(rs.values())  # number of roots

    # get taylor coefficients
    f = sympy.utilities.lambdify(list(g_z.free_symbols), g_z)
    taylor_coeffs = mpmath.taylor(f, 0, d - 1)  # get the first d terms of taylor expansion

    #
    # calculate path_n
    #

    n = sympy.var('n')
    e_path_n = 0
    e_upper_n = 0
    coeff = []

    for i in xrange(1, D + 1):
        ri, mi = rs.items()[i - 1]
        for j in xrange(mi):
            c_ij = sympy.var('c_' + str(i) + str(j))
            coeff.append(c_ij)
            e_path_n += c_ij * (n ** j) * ((1 / ri) ** n)
            if ri.is_complex:
                ri = sympy.functions.Abs(ri)
            e_upper_n += c_ij * (n ** j) * ((1 / ri) ** n)
    equations = []

    for i, c in enumerate(taylor_coeffs):
        equations.append(sympy.Eq(e_path_n.subs(n, i), c))
    coeff_sol = sympy.linsolve(equations, coeff)

    # assert unique solution
    assert type(coeff_sol) == sympy.sets.FiniteSet, "Zero or more solutions returned for path_n coefficients"
    coeff_sol = list(coeff_sol)[0]
    coeff_sol = [sympy.N(c, ROUND) for c in coeff_sol]

    for val, var in zip(coeff_sol, coeff):
        name = var.name
        e_path_n = e_path_n.subs(name, val)
        e_upper_n = e_upper_n.subs(name, val)

    return sympy.utilities.lambdify(list(e_path_n.free_symbols), e_path_n), sympy.utilities.lambdify(
        list(e_upper_n.free_symbols), e_upper_n)