Python sympy 模块，DeferredVector() 实例源码

def _integrate_exact(f, quadrilateral):
    xi = sympy.DeferredVector('xi')
    pxi = quadrilateral[0] * 0.25*(1.0 + xi[0])*(1.0 + xi[1]) \
        + quadrilateral[1] * 0.25*(1.0 - xi[0])*(1.0 + xi[1]) \
        + quadrilateral[2] * 0.25*(1.0 - xi[0])*(1.0 - xi[1]) \
        + quadrilateral[3] * 0.25*(1.0 + xi[0])*(1.0 - xi[1])
    pxi = [
        sympy.expand(pxi[0]),
        sympy.expand(pxi[1]),
        ]
    # determinant of the transformation matrix
    det_J = \
        + sympy.diff(pxi[0], xi[0]) * sympy.diff(pxi[1], xi[1]) \
        - sympy.diff(pxi[1], xi[0]) * sympy.diff(pxi[0], xi[1])
    # we cannot use abs(), see <https://github.com/sympy/sympy/issues/4212>.
    abs_det_J = sympy.Piecewise((det_J, det_J >= 0), (-det_J, det_J < 0))

    g_xi = f(pxi)

    exact = sympy.integrate(
        sympy.integrate(abs_det_J * g_xi, (xi[1], -1, 1)),
        (xi[0], -1, 1)
        )
    return float(exact)

def _integrate_exact(f, tetrahedron):
    #
    # Note that
    #
    #     \int_T f(x) dx = \int_T0 |J(xi)| f(P(xi)) dxi
    #
    # with
    #
    #     P(xi) = x0 * (1-xi[0]-xi[1]) + x1 * xi[0] + x2 * xi[1].
    #
    # and T0 being the reference tetrahedron [(0.0, 0.0), (1.0, 0.0), (0.0,
    # 1.0)].
    # The determinant of the transformation matrix J equals twice the volume of
    # the tetrahedron. (See, e.g.,
    # <http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF>).
    #
    xi = sympy.DeferredVector('xi')
    x_xi = \
        + tetrahedron[0] * (1 - xi[0] - xi[1] - xi[2]) \
        + tetrahedron[1] * xi[0] \
        + tetrahedron[2] * xi[1] \
        + tetrahedron[3] * xi[2]
    abs_det_J = 6 * quadpy.tetrahedron.volume(tetrahedron)
    exact = sympy.integrate(
        sympy.integrate(
          sympy.integrate(abs_det_J * f(x_xi), (xi[2], 0, 1-xi[0]-xi[1])),
          (xi[1], 0, 1-xi[0])
          ),
        (xi[0], 0, 1)
      )
    return float(exact)

def _integrate_exact(f, triangle):
    #
    # Note that
    #
    #     \int_T f(x) dx = \int_T0 |J(xi)| f(P(xi)) dxi
    #
    # with
    #
    #     P(xi) = x0 * (1-xi[0]-xi[1]) + x1 * xi[0] + x2 * xi[1].
    #
    # and T0 being the reference triangle [(0.0, 0.0), (1.0, 0.0), (0.0,
    # 1.0)].
    # The determinant of the transformation matrix J equals twice the volume of
    # the triangle. (See, e.g.,
    # <http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF>).
    #
    xi = sympy.DeferredVector('xi')
    x_xi = \
        + triangle[0] * (1 - xi[0] - xi[1]) \
        + triangle[1] * xi[0] \
        + triangle[2] * xi[1]
    abs_det_J = 2 * quadpy.triangle.volume(triangle)
    exact = sympy.integrate(
        sympy.integrate(abs_det_J * f(x_xi), (xi[1], 0, 1-xi[0])),
        (xi[0], 0, 1)
        )
    return float(exact)

def test_lambdify_matrix_vec_input():
    X = sympy.DeferredVector('X')
    M = Matrix([
        [X[0]**2, X[0]*X[1], X[0]*X[2]],
        [X[1]*X[0], X[1]**2, X[1]*X[2]],
        [X[2]*X[0], X[2]*X[1], X[2]**2]])
    f = lambdify(X, M, "numpy")

    Xh = array([1.0, 2.0, 3.0])
    expected = matrix([[Xh[0]**2, Xh[0]*Xh[1], Xh[0]*Xh[2]],
                       [Xh[1]*Xh[0], Xh[1]**2, Xh[1]*Xh[2]],
                       [Xh[2]*Xh[0], Xh[2]*Xh[1], Xh[2]**2]])
    actual = f(Xh)
    assert numpy.allclose(actual, expected)

def test_lambdify_matrix_vec_input():
    X = sympy.DeferredVector('X')
    M = Matrix([
        [X[0]**2, X[0]*X[1], X[0]*X[2]],
        [X[1]*X[0], X[1]**2, X[1]*X[2]],
        [X[2]*X[0], X[2]*X[1], X[2]**2]])
    f = lambdify(X, M, [{'ImmutableMatrix': numpy.array}, "numpy"])

    Xh = array([1.0, 2.0, 3.0])
    expected = array([[Xh[0]**2, Xh[0]*Xh[1], Xh[0]*Xh[2]],
                      [Xh[1]*Xh[0], Xh[1]**2, Xh[1]*Xh[2]],
                      [Xh[2]*Xh[0], Xh[2]*Xh[1], Xh[2]**2]])
    actual = f(Xh)
    assert numpy.allclose(actual, expected)

def _integrate_exact(k, pyra):
    def f(x):
        return x[0]**int(k[0]) * x[1]**int(k[1]) * x[2]**int(k[2])

    # map the reference hexahedron [-1,1]^3 to the pyramid
    xi = sympy.DeferredVector('xi')
    pxi = (
        + pyra[0] * (1 - xi[0])*(1 - xi[1])*(1 - xi[2]) / 8
        + pyra[1] * (1 + xi[0])*(1 - xi[1])*(1 - xi[2]) / 8
        + pyra[2] * (1 + xi[0])*(1 + xi[1])*(1 - xi[2]) / 8
        + pyra[3] * (1 - xi[0])*(1 + xi[1])*(1 - xi[2]) / 8
        + pyra[4] * (1 + xi[2]) / 2
        )

    pxi = [
        sympy.expand(pxi[0]),
        sympy.expand(pxi[1]),
        sympy.expand(pxi[2]),
        ]
    # determinant of the transformation matrix
    J = sympy.Matrix([
        [sympy.diff(pxi[0], xi[0]),
         sympy.diff(pxi[0], xi[1]),
         sympy.diff(pxi[0], xi[2])],
        [sympy.diff(pxi[1], xi[0]),
         sympy.diff(pxi[1], xi[1]),
         sympy.diff(pxi[1], xi[2])],
        [sympy.diff(pxi[2], xi[0]),
         sympy.diff(pxi[2], xi[1]),
         sympy.diff(pxi[2], xi[2])],
        ])
    det_J = sympy.det(J)
    # we cannot use abs(), see <https://github.com/sympy/sympy/issues/4212>.
    # abs_det_J = sympy.Piecewise((det_J, det_J >= 0), (-det_J, det_J < 0))
    # This is quite the leap of faith, but sympy will cowardly bail out
    # otherwise.
    abs_det_J = det_J

    exact = sympy.integrate(
        sympy.integrate(
            sympy.integrate(abs_det_J * f(pxi), (xi[2], -1, 1)),
            (xi[1], -1, +1)
            ),
        (xi[0], -1, +1)
        )

    return float(exact)

def _integrate_exact(f, hexa):
    xi = sympy.DeferredVector('xi')
    pxi = \
        + hexa[0] * 0.125*(1.0 - xi[0])*(1.0 - xi[1])*(1.0 - xi[2]) \
        + hexa[1] * 0.125*(1.0 + xi[0])*(1.0 - xi[1])*(1.0 - xi[2]) \
        + hexa[2] * 0.125*(1.0 + xi[0])*(1.0 + xi[1])*(1.0 - xi[2]) \
        + hexa[3] * 0.125*(1.0 - xi[0])*(1.0 + xi[1])*(1.0 - xi[2]) \
        + hexa[4] * 0.125*(1.0 - xi[0])*(1.0 - xi[1])*(1.0 + xi[2]) \
        + hexa[5] * 0.125*(1.0 + xi[0])*(1.0 - xi[1])*(1.0 + xi[2]) \
        + hexa[6] * 0.125*(1.0 + xi[0])*(1.0 + xi[1])*(1.0 + xi[2]) \
        + hexa[7] * 0.125*(1.0 - xi[0])*(1.0 + xi[1])*(1.0 + xi[2])
    pxi = [
        sympy.expand(pxi[0]),
        sympy.expand(pxi[1]),
        sympy.expand(pxi[2]),
        ]
    # determinant of the transformation matrix
    J = sympy.Matrix([
        [sympy.diff(pxi[0], xi[0]),
         sympy.diff(pxi[0], xi[1]),
         sympy.diff(pxi[0], xi[2])],
        [sympy.diff(pxi[1], xi[0]),
         sympy.diff(pxi[1], xi[1]),
         sympy.diff(pxi[1], xi[2])],
        [sympy.diff(pxi[2], xi[0]),
         sympy.diff(pxi[2], xi[1]),
         sympy.diff(pxi[2], xi[2])],
        ])
    det_J = sympy.det(J)
    # we cannot use abs(), see <https://github.com/sympy/sympy/issues/4212>.
    abs_det_J = sympy.Piecewise((det_J, det_J >= 0), (-det_J, det_J < 0))
    g_xi = f(pxi)
    exact = \
        sympy.integrate(
            sympy.integrate(
                sympy.integrate(abs_det_J * g_xi, (xi[2], -1, 1)),
                (xi[1], -1, 1)
            ),
            (xi[0], -1, 1)
        )
    return float(exact)


# pylint: disable=too-many-arguments

def _newton_cotes(n, point_fun):
    '''
    Construction after

    P. Silvester,
    Symmetric quadrature formulae for simplexes
    Math. Comp., 24, 95-100 (1970),
    <https://doi.org/10.1090/S0025-5718-1970-0258283-6>.
    '''
    degree = n

    # points
    idx = numpy.array([
        [i, j, n-i-j]
        for i in range(n + 1)
        for j in range(n + 1 - i)
        ])
    bary = point_fun(idx, n)

    # weights
    if n == 0:
        weights = numpy.ones(1)
        return bary, weights, degree

    def get_poly(t, m, n):
        return sympy.prod([
            sympy.poly(
                (t - point_fun(k, n)) / (point_fun(m, n) - point_fun(k, n))
                )
            for k in range(m)
            ])
    weights = numpy.empty(len(bary))
    idx = 0
    for i in range(n + 1):
        for j in range(n + 1 - i):
            k = n - i - j
            # Define the polynomial which to integrate over the
            # tetrahedron.
            t = sympy.DeferredVector('t')
            g = get_poly(t[0], i, n) \
                * get_poly(t[1], j, n) \
                * get_poly(t[2], k, n)
            # The integral of monomials over a tetrahedron are well-known,
            # see Silvester.
            weights[idx] = numpy.sum([
                 c * numpy.prod([math.factorial(l) for l in m]) * 2
                 / math.factorial(numpy.sum(m) + 2)
                 for m, c in zip(g.monoms(), g.coeffs())
                 ])
            idx += 1
    return bary, weights, degree

def __init__(self,
                f,
                past,
                helpers = (),
                control_pars = (),
                n_basic = None,
                tangent_indices = (),
            ):
        self.past = past
        self.t, self.y, self.diff = self.past[-1]
        self.n = len(self.y)
        self.n_basic = n_basic or self.n
        self.tangent_indices = tangent_indices
        self.last_garbage = -1
        self.old_new_y = None

        self.parameters = []

        from jitcdde._jitcdde import t, y, past_y, anchors
        from sympy import DeferredVector, sympify, lambdify
        Y = DeferredVector("Y")
        substitutions = list(reversed(helpers)) + [(y(i),Y[i]) for i in range(self.n)]

        past_calls = 0
        f_wc = []
        for entry in f():
            new_entry = sympify(entry).subs(substitutions).simplify(ratio=1.0)
            past_calls += new_entry.count(anchors)
            f_wc.append(new_entry)

        F = lambdify(
            [t, Y] + list(control_pars),
            f_wc,
            [
                {
                    anchors.name: self.get_past_anchors,
                    past_y .name: interpolate
                },
                "math"
            ]
            )

        self.f = lambda *args: np.array(F(*args)).flatten()

        self.anchor_mem = (len(past)-1)*np.ones(past_calls, dtype=int)