Python sympy.integrate() Examples

def test_exp():
    """
    Make sure that symfit.distributions.Exp produces the expected
    sympy expression.
    """
    l = Parameter(positive=True)
    x = Variable()

    new = l * sympy.exp(- l * x)
    assert isinstance(new, sympy.Expr)
    e = Exp(x, l)
    assert issubclass(e.__class__, sympy.Expr)
    assert new == e

    # A pdf should always integrate to 1 on its domain
    assert sympy.integrate(e, (x, 0, sympy.oo)) == 1 

def test_gaussian():
    """
    Make sure that symfit.distributions.Gaussians produces the expected
    sympy expression.
    """
    x0 = Parameter()
    sig = Parameter(positive=True)
    x = Variable()

    new = sympy.exp(-(x - x0)**2/(2*sig**2))/sympy.sqrt((2*sympy.pi*sig**2))
    assert isinstance(new, sympy.Expr)
    g = Gaussian(x, x0, sig)
    assert issubclass(g.__class__, sympy.Expr)
    assert new == g

    # A pdf should always integrate to 1 on its domain
    assert sympy.integrate(g, (x, -sympy.oo, sympy.oo)) == 1 

def test_scheme(scheme):
    assert scheme.points.dtype in [numpy.float64, numpy.int64], scheme.name
    assert scheme.weights.dtype in [numpy.float64, numpy.int64], scheme.name

    print(scheme)

    # Test integration until we get to a polynomial degree `d` that can no longer be
    # integrated exactly. The scheme's degree is `d-1`.
    pyra = numpy.array(
        [[-1, -1, -1], [+1, -1, -1], [+1, +1, -1], [-1, +1, -1], [0, 0, 1]]
    )
    degree, err = check_degree(
        lambda poly: scheme.integrate(poly, pyra),
        lambda k: _integrate_exact(k, pyra),
        3,
        scheme.degree + 1,
        tol=scheme.test_tolerance,
    )
    assert (
        degree >= scheme.degree
    ), "{} -- Observed: {}, expected: {} (max err: {:.3e})".format(
        scheme.name, degree, scheme.degree, err
    ) 

def test_scheme(scheme):
    assert scheme.points.dtype in [numpy.float64, numpy.int64], scheme.name
    assert scheme.weights.dtype in [numpy.float64, numpy.int64], scheme.name

    print(scheme)

    # Test integration until we get to a polynomial degree `d` that can no
    # longer be integrated exactly. The scheme's degree is `d-1`.
    t3 = numpy.array(
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
    )

    degree, err = check_degree(
        lambda poly: scheme.integrate(poly, t3),
        integrate_monomial_over_unit_simplex,
        3,
        scheme.degree + 1,
        scheme.test_tolerance,
    )

    assert (
        degree >= scheme.degree
    ), "{} -- observed: {}, expected: {} (max err: {:.3e})".format(
        scheme.name, degree, scheme.degree, err
    ) 

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.t2.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 newton_cotes_open(index, **kwargs):
    """
    Open Newton-Cotes formulae.
    <https://math.stackexchange.com/a/1959071/36678>
    """
    points = numpy.linspace(-1.0, 1.0, index + 2)[1:-1]
    degree = index if (index + 1) % 2 == 0 else index - 1
    #
    n = index + 1
    weights = numpy.empty(n - 1)
    t = sympy.Symbol("t")
    for r in range(1, n):
        # Compare with get_weights().
        f = sympy.prod([(t - i) for i in range(1, n) if i != r])
        alpha = (
            2
            * (-1) ** (n - r + 1)
            * sympy.integrate(f, (t, 0, n), **kwargs)
            / (math.factorial(r - 1) * math.factorial(n - 1 - r))
            / n
        )
        weights[r - 1] = alpha
    return C1Scheme("Newton-Cotes (open)", degree, weights, points) 

def integrate(f, a, b, **kwargs):
    """Symbolically calculate the integrals

      int_a^b f_k(x) dx.

    Useful for computing the moments `w(x) * P_k(x)`, e.g.,

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

    Any keyword arguments are passed directly to `sympy.integrate` function.
    """
    x = sympy.Symbol("x")
    return numpy.array([sympy.integrate(fun, (x, a, b), **kwargs) for fun in f(x)]) 

def stieltjes(w, a, b, n, **kwargs):
    t = sympy.Symbol("t")

    alpha = n * [None]
    beta = n * [None]
    mu = n * [None]
    pi = n * [None]

    k = 0
    pi[k] = 1
    mu[k] = sympy.integrate(pi[k] ** 2 * w(t), (t, a, b), **kwargs)
    alpha[k] = sympy.integrate(t * pi[k] ** 2 * w(t), (t, a, b), **kwargs) / mu[k]
    beta[k] = mu[0]  # not used, by convention mu[0]

    k = 1
    pi[k] = (t - alpha[k - 1]) * pi[k - 1]
    mu[k] = sympy.integrate(pi[k] ** 2 * w(t), (t, a, b), **kwargs)
    alpha[k] = sympy.integrate(t * pi[k] ** 2 * w(t), (t, a, b), **kwargs) / mu[k]
    beta[k] = mu[k] / mu[k - 1]

    for k in range(2, n):
        pi[k] = (t - alpha[k - 1]) * pi[k - 1] - beta[k - 1] * pi[k - 2]
        mu[k] = sympy.integrate(pi[k] ** 2 * w(t), (t, a, b), **kwargs)
        alpha[k] = sympy.integrate(t * pi[k] ** 2 * w(t), (t, a, b), **kwargs) / mu[k]
        beta[k] = mu[k] / mu[k - 1]

    return alpha, beta 

def _integrate_exact(f, t3):
    #
    # 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 t3 [(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 t3. (See, e.g.,
    # <http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF>).
    #
    xi = sympy.DeferredVector("xi")
    x_xi = (
        +t3[0] * (1 - xi[0] - xi[1] - xi[2])
        + t3[1] * xi[0]
        + t3[2] * xi[1]
        + t3[3] * xi[2]
    )
    abs_det_J = 6 * quadpy.t3.volume(t3)
    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 test_scheme(scheme):
    # Test integration until we get to a polynomial degree `d` that can no longer be
    # integrated exactly. The scheme's degree is `d-1`.
    assert scheme.points.dtype in [numpy.float64, numpy.int64], scheme.name
    assert scheme.weights.dtype in [numpy.float64, numpy.int64], scheme.name

    print(scheme)

    def eval_orthopolys(x):
        return numpy.concatenate(
            orthopy.quadrilateral.tree(x, scheme.degree + 1, symbolic=False)
        )

    quad = quadpy.c2.rectangle_points([-1.0, +1.0], [-1.0, +1.0])
    vals = scheme.integrate(eval_orthopolys, quad)
    # Put vals back into the tree structure:
    # len(approximate[k]) == k+1
    approximate = [
        vals[k * (k + 1) // 2 : (k + 1) * (k + 2) // 2]
        for k in range(scheme.degree + 2)
    ]

    exact = [numpy.zeros(k + 1) for k in range(scheme.degree + 2)]
    exact[0][0] = 2.0

    degree, err = check_degree_ortho(approximate, exact, abs_tol=scheme.test_tolerance)

    assert (
        degree >= scheme.degree
    ), "{} -- observed: {}, expected: {} (max err: {:.3e})".format(
        scheme.name, degree, scheme.degree, err
    ) 

def _integrate_exact(f, c2):
    xi = sympy.DeferredVector("xi")
    pxi = (
        c2[0] * 0.25 * (1.0 + xi[0]) * (1.0 + xi[1])
        + c2[1] * 0.25 * (1.0 - xi[0]) * (1.0 + xi[1])
        + c2[2] * 0.25 * (1.0 - xi[0]) * (1.0 - xi[1])
        + c2[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 phi(i, t):
    if i == 0:
        return 1  # Initial history x(t)=1 on [-1,0]
    else:
        return phi(i-1, i-1) - integrate(phi(i-1, xi-1), (xi, i-1, t)) 

def newton_cotes_closed(index, **kwargs):
    """
    Closed Newton-Cotes formulae.
    <https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas#Closed_Newton.E2.80.93Cotes_formulae>,
    <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
    """
    points = numpy.linspace(-1.0, 1.0, index + 1)
    degree = index + 1 if index % 2 == 0 else index

    # Formula (26) from
    # <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
    # Note that Sympy carries out all operations in rationals, i.e.,
    # _exactly_. Only at the end, the rational is converted into a float.
    n = index
    weights = numpy.empty(n + 1)
    t = sympy.Symbol("t")
    for r in range(n + 1):
        # Compare with get_weights().
        f = sympy.prod([(t - i) for i in range(n + 1) if i != r])
        alpha = (
            2
            * (-1) ** (n - r)
            * sympy.integrate(f, (t, 0, n), **kwargs)
            / (math.factorial(r) * math.factorial(n - r))
            / index
        )
        weights[r] = alpha
    return C1Scheme("Newton-Cotes (closed)", degree, weights, points) 

def quadrature_degree(self, u, h, **kwargs):
        r"""Return the quadrature degree necessary to integrate the action
        functional accurately

        Firedrake uses a very conservative algorithm for estimating the
        number of quadrature points necessary to integrate a given
        expression. By exploiting known structure of the problem, we can
        reduce the number of quadrature points while preserving accuracy.
        """
        xdegree_u, zdegree_u = u.ufl_element().degree()
        degree_h = h.ufl_element().degree()[0]
        return (3 * (xdegree_u - 1) + 2 * degree_h,
                3 * max(zdegree_u - 1, 0) + zdegree_u + 1) 

def _pressure_approx(N):
    ζ, ζ_sl = sympy.symbols('ζ ζ_sl', real=True, positive=True)

    def coefficient(n):
        Sn = _legendre(n, ζ)
        norm_square = sympy.integrate(Sn**2, (ζ, 0, 1))
        return sympy.integrate((ζ_sl - ζ) * Sn, (ζ, 0, ζ_sl)) / norm_square

    polynomial = sum([coefficient(n) * _legendre(n, ζ) for n in range(N)])
    return sympy.lambdify((ζ, ζ_sl), sympy.simplify(polynomial)) 

def sol(self):
        # Do the inner product! - we are in cyl coordinates!
        j, Sig = self.fcts()
        jTSj = j.T*Sig*j
        # we are integrating in cyl coordinates
        ans = sympy.integrate(
            sympy.integrate(
                sympy.integrate(r * jTSj, (r, 0, 1)),
                (t, 0, 2*sympy.pi)
            ),
            (z, 0, 1)
        )[0] # The `[0]` is to make it a number rather than a matrix

        return ans 

def sol(self):
        h, Sig = self.fcts()

        hTSh = h.T*Sig*h
        ans  = sympy.integrate(
            sympy.integrate(
                sympy.integrate(r * hTSh, (r, 0, 1)),
                (t, 0, 2*sympy.pi)),
            (z, 0, 1)
        )[0] # The `[0]` is to make it a scalar

        return ans 

def sol(self):
        # Do the inner product! - we are in cyl coordinates!
        j, Sig = self.fcts()
        jTSj = j.T*Sig*j
        # we are integrating in cyl coordinates
        ans  = sympy.integrate(sympy.integrate(sympy.integrate(r * jTSj,
                                                               (r, 0, 1)),
                                               (t, 0, 2*sympy.pi)),
                                (z, 0, 1))[0] # The `[0]` is to make it an int.

        return ans 

def sol(self):
        h, Sig = self.fcts()
        # Do the inner product! - we are in cyl coordinates!
        hTSh = h.T*Sig*h
        ans  = sympy.integrate(sympy.integrate(sympy.integrate(r * hTSh,
                                                               (r, 0, 1)),
                                               (t, 0, 2*sympy.pi)),
                                (z, 0, 1))[0] # The `[0]` is to make it an int.
        return ans 

def integral(self) -> Dict[ChannelID, ExpressionScalar]:
        return {self.__channel: ExpressionScalar(
            sympy.integrate(self.__expression.sympified_expression, ('t', 0, self.duration.sympified_expression))
        )} 

def test_multidim():
    scheme = quadpy.t2.dunavant_05()

    numpy.random.seed(0)
    # simple scalar integration
    tri = numpy.random.rand(3, 2)
    val = scheme.integrate(lambda x: numpy.sin(x[0]), tri)
    assert val.shape == ()

    # scalar integration on 4 subdomains
    tri = numpy.random.rand(3, 4, 2)
    val = scheme.integrate(lambda x: numpy.sin(x[0]), tri)
    assert val.shape == (4,)

    # scalar integration in 4D
    tri = numpy.random.rand(3, 4)
    val = scheme.integrate(lambda x: numpy.sin(x[0]), tri)
    assert val.shape == ()

    # vector-valued integration on 4 subdomains
    tri = numpy.random.rand(3, 4, 2)
    val = scheme.integrate(lambda x: [numpy.sin(x[0]), numpy.cos(x[1])], tri)
    assert val.shape == (2, 4)

    # vector-valued integration in 4D
    tri = numpy.random.rand(3, 4)
    val = scheme.integrate(lambda x: [numpy.sin(x[0]), numpy.cos(x[1])], tri)
    assert val.shape == (2,)

    # # another vector-valued integration in 3D
    # # This is one case where the integration routine may not properly recognize the
    # # dimensionality of the domain. Use the `dim` parameter.
    # val = scheme.integrate(
    #     lambda x: [
    #         x[0] + numpy.sin(x[1]),
    #         numpy.cos(x[0]) * x[2],
    #         numpy.sin(x[0]) + x[1] + x[2],
    #     ],
    #     [[0.0, 1.0, 2.0], [1.0, 2.0, 3.0]],
    #     dim=1,
    # )
    # assert val.shape == (3,) 

def test_scheme(scheme, print_degree=False):
    assert scheme.points.dtype in [numpy.float64, numpy.int64], scheme.name
    assert scheme.weights.dtype in [numpy.float64, numpy.int64], scheme.name

    print(scheme)

    x = [-1.0, +1.0]
    y = [-1.0, +1.0]
    z = [-1.0, +1.0]
    hexa = quadpy.c3.cube_points(x, y, z)

    # degree = check_degree(
    #     lambda poly: scheme.integrate(poly, hexa),
    #     lambda k: _integrate_exact2(k, x[0], x[1], y[0], y[1], z[0], z[1]),
    #     3,
    #     scheme.degree + 1,
    #     tol=tol,
    # )
    # if print_degree:
    #     print("Detected degree {}, scheme degree {}.".format(degree, scheme.degree))
    # assert degree == scheme.degree, scheme.name

    def eval_orthopolys(x):
        return numpy.concatenate(
            orthopy.hexahedron.tree(x, scheme.degree + 1, symbolic=False)
        )

    vals = scheme.integrate(eval_orthopolys, hexa)
    # Put vals back into the tree structure:
    # len(approximate[k]) == k+1
    approximate = [
        vals[k * (k + 1) * (k + 2) // 6 : (k + 1) * (k + 2) * (k + 3) // 6]
        for k in range(scheme.degree + 2)
    ]

    exact = [numpy.zeros(len(s)) for s in approximate]
    exact[0][0] = numpy.sqrt(2.0) * 2

    degree, err = check_degree_ortho(approximate, exact, abs_tol=scheme.test_tolerance)

    assert (
        degree >= scheme.degree
    ), "{} -- Observed: {}, expected: {} (max err: {:.3e})".format(
        scheme.name, degree, scheme.degree, err
    ) 

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 p3
    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 generate_calculus_integrate_sample(vlist, ops, min_depth, max_depth,
                                       functions):
  """Randomly generate a symbolic integral dataset sample.

  Given an input expression, produce the indefinite integral.

  See go/symbolic-math-dataset.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.
    functions: Defines special functions. A dict mapping human readable string
        names, like "log", "exp", "sin", "cos", etc., to single chars. Each
        function gets a unique token, like "L" for "log".

  Returns:
    sample: String representation of the input. Will be of the form
        'var:expression'.
    target: String representation of the solution.
  """
  var_index = random.randrange(len(vlist))
  var = vlist[var_index]
  consts = vlist[:var_index] + vlist[var_index + 1:]

  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr_with_required_var(depth, var, consts, ops)

  expr_str = str(expr)
  sample = var + ":" + expr_str
  target = format_sympy_expr(
      sympy.integrate(expr_str, sympy.Symbol(var)), functions=functions)
  return sample, target


# AlgebraConfig holds objects required to generate the algebra inverse
# dataset. See go/symbolic-math-dataset.
# vlist: Variable list. A list of chars.
# dlist: Numberical digit list. A list of chars.
# flist: List of special function names. A list of chars.
# functions: Dict of special function names. Maps human readable string names to
#     single char names used in flist.
# ops: Dict mapping op symbols (chars) to ExprOp instances.
# solve_ops: Encodes rules for how to algebraically cancel out each operation.
#     See doc-string for `algebra_inverse_solve`.
# int_encoder: Function that maps a string to a list of tokens. Use this to
#     encode an expression to feed into a model.
# int_decoder: Function that maps a list of tokens to a string. Use this to
#     convert model input or output into a human readable string. 

def calculus_integrate(alphabet_size=26,
                       min_depth=0,
                       max_depth=2,
                       nbr_cases=10000):
  """Generate the calculus integrate dataset.

  Each sample is a symbolic math expression involving unknown variables. The
  task is to take the indefinite integral of the expression. The target is the
  resulting expression.

  Args:
    alphabet_size: How many possible variables there are. Max 26.
    min_depth: Minimum depth of the expression trees on both sides of the
        equals sign in the equation.
    max_depth: Maximum depth of the expression trees on both sides of the
        equals sign in the equation.
    nbr_cases: The number of cases to generate.

  Yields:
    A dictionary {"inputs": input-list, "targets": target-list} where
    input-list are the tokens encoding the variable to integrate with respect
    to and the expression to integrate, and target-list is a list of tokens
    encoding the resulting math expression after integrating.

  Raises:
    ValueError: If `max_depth` < `min_depth`, or if alphabet_size > 26.
  """
  if max_depth < min_depth:
    raise ValueError("max_depth must be greater than or equal to min_depth. "
                     "Got max_depth=%s, min_depth=%s" % (max_depth, min_depth))

  # Don't allow alphabet to use capital letters. Those are reserved for function
  # names.
  if alphabet_size > 26:
    raise ValueError(
        "alphabet_size must not be greater than 26. Got %s." % alphabet_size)

  functions = {"log": "L"}
  alg_cfg = math_dataset_init(alphabet_size, digits=5, functions=functions)
  nbr_case = 0
  while nbr_case < nbr_cases:
    try:
      sample, target = generate_calculus_integrate_sample(
          alg_cfg.vlist,
          list(alg_cfg.ops.values()), min_depth, max_depth, alg_cfg.functions)
      yield {
          "inputs": alg_cfg.int_encoder(sample),
          "targets": alg_cfg.int_encoder(target)
      }
    except:  # pylint:disable=bare-except
      continue
    if nbr_case % 10000 == 0:
      print(" calculus_integrate: generating case %d." % nbr_case)
    nbr_case += 1 

def generate_calculus_integrate_sample(vlist, ops, min_depth, max_depth,
                                       functions):
  """Randomly generate a symbolic integral dataset sample.

  Given an input expression, produce the indefinite integral.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.
    functions: Defines special functions. A dict mapping human readable string
        names, like "log", "exp", "sin", "cos", etc., to single chars. Each
        function gets a unique token, like "L" for "log".

  Returns:
    sample: String representation of the input. Will be of the form
        'var:expression'.
    target: String representation of the solution.
  """
  var_index = random.randrange(len(vlist))
  var = vlist[var_index]
  consts = vlist[:var_index] + vlist[var_index + 1:]

  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr_with_required_var(depth, var, consts, ops)

  expr_str = str(expr)
  sample = var + ":" + expr_str
  target = format_sympy_expr(
      sympy.integrate(expr_str, sympy.Symbol(var)), functions=functions)
  return sample, target


# AlgebraConfig holds objects required to generate the algebra inverse
# dataset.
# vlist: Variable list. A list of chars.
# dlist: Numberical digit list. A list of chars.
# flist: List of special function names. A list of chars.
# functions: Dict of special function names. Maps human readable string names to
#     single char names used in flist.
# ops: Dict mapping op symbols (chars) to ExprOp instances.
# solve_ops: Encodes rules for how to algebraically cancel out each operation.
#     See doc-string for `algebra_inverse_solve`.
# int_encoder: Function that maps a string to a list of tokens. Use this to
#     encode an expression to feed into a model.
# int_decoder: Function that maps a list of tokens to a string. Use this to
#     convert model input or output into a human readable string. 

def calculus_integrate(alphabet_size=26,
                       min_depth=0,
                       max_depth=2,
                       nbr_cases=10000):
  """Generate the calculus integrate dataset.

  Each sample is a symbolic math expression involving unknown variables. The
  task is to take the indefinite integral of the expression. The target is the
  resulting expression.

  Args:
    alphabet_size: How many possible variables there are. Max 26.
    min_depth: Minimum depth of the expression trees on both sides of the
        equals sign in the equation.
    max_depth: Maximum depth of the expression trees on both sides of the
        equals sign in the equation.
    nbr_cases: The number of cases to generate.

  Yields:
    A dictionary {"inputs": input-list, "targets": target-list} where
    input-list are the tokens encoding the variable to integrate with respect
    to and the expression to integrate, and target-list is a list of tokens
    encoding the resulting math expression after integrating.

  Raises:
    ValueError: If `max_depth` < `min_depth`, or if alphabet_size > 26.
  """
  if max_depth < min_depth:
    raise ValueError("max_depth must be greater than or equal to min_depth. "
                     "Got max_depth=%s, min_depth=%s" % (max_depth, min_depth))

  # Don't allow alphabet to use capital letters. Those are reserved for function
  # names.
  if alphabet_size > 26:
    raise ValueError(
        "alphabet_size must not be greater than 26. Got %s." % alphabet_size)

  functions = {"log": "L"}
  alg_cfg = math_dataset_init(alphabet_size, digits=5, functions=functions)
  nbr_case = 0
  while nbr_case < nbr_cases:
    try:
      sample, target = generate_calculus_integrate_sample(
          alg_cfg.vlist,
          list(alg_cfg.ops.values()), min_depth, max_depth, alg_cfg.functions)
      yield {
          "inputs": alg_cfg.int_encoder(sample),
          "targets": alg_cfg.int_encoder(target)
      }
    except:  # pylint:disable=bare-except
      continue
    if nbr_case % 10000 == 0:
      print(" calculus_integrate: generating case %d." % nbr_case)
    nbr_case += 1 

def generate_calculus_integrate_sample(vlist, ops, min_depth, max_depth,
                                       functions):
  """Randomly generate a symbolic integral dataset sample.

  Given an input expression, produce the indefinite integral.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.
    functions: Defines special functions. A dict mapping human readable string
        names, like "log", "exp", "sin", "cos", etc., to single chars. Each
        function gets a unique token, like "L" for "log".

  Returns:
    sample: String representation of the input. Will be of the form
        'var:expression'.
    target: String representation of the solution.
  """
  var_index = random.randrange(len(vlist))
  var = vlist[var_index]
  consts = vlist[:var_index] + vlist[var_index + 1:]

  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr_with_required_var(depth, var, consts, ops)

  expr_str = str(expr)
  sample = var + ":" + expr_str
  target = format_sympy_expr(
      sympy.integrate(expr_str, sympy.Symbol(var)), functions=functions)
  return sample, target


# AlgebraConfig holds objects required to generate the algebra inverse
# dataset.
# vlist: Variable list. A list of chars.
# dlist: Numberical digit list. A list of chars.
# flist: List of special function names. A list of chars.
# functions: Dict of special function names. Maps human readable string names to
#     single char names used in flist.
# ops: Dict mapping op symbols (chars) to ExprOp instances.
# solve_ops: Encodes rules for how to algebraically cancel out each operation.
#     See doc-string for `algebra_inverse_solve`.
# int_encoder: Function that maps a string to a list of tokens. Use this to
#     encode an expression to feed into a model.
# int_decoder: Function that maps a list of tokens to a string. Use this to
#     convert model input or output into a human readable string. 

def calculus_integrate(alphabet_size=26,
                       min_depth=0,
                       max_depth=2,
                       nbr_cases=10000):
  """Generate the calculus integrate dataset.

  Each sample is a symbolic math expression involving unknown variables. The
  task is to take the indefinite integral of the expression. The target is the
  resulting expression.

  Args:
    alphabet_size: How many possible variables there are. Max 26.
    min_depth: Minimum depth of the expression trees on both sides of the
        equals sign in the equation.
    max_depth: Maximum depth of the expression trees on both sides of the
        equals sign in the equation.
    nbr_cases: The number of cases to generate.

  Yields:
    A dictionary {"inputs": input-list, "targets": target-list} where
    input-list are the tokens encoding the variable to integrate with respect
    to and the expression to integrate, and target-list is a list of tokens
    encoding the resulting math expression after integrating.

  Raises:
    ValueError: If `max_depth` < `min_depth`, or if alphabet_size > 26.
  """
  if max_depth < min_depth:
    raise ValueError("max_depth must be greater than or equal to min_depth. "
                     "Got max_depth=%s, min_depth=%s" % (max_depth, min_depth))

  # Don't allow alphabet to use capital letters. Those are reserved for function
  # names.
  if alphabet_size > 26:
    raise ValueError(
        "alphabet_size must not be greater than 26. Got %s." % alphabet_size)

  functions = {"log": "L"}
  alg_cfg = math_dataset_init(alphabet_size, digits=5, functions=functions)
  nbr_case = 0
  while nbr_case < nbr_cases:
    try:
      sample, target = generate_calculus_integrate_sample(
          alg_cfg.vlist,
          list(alg_cfg.ops.values()), min_depth, max_depth, alg_cfg.functions)
      yield {
          "inputs": alg_cfg.int_encoder(sample),
          "targets": alg_cfg.int_encoder(target)
      }
    except:  # pylint:disable=bare-except
      continue
    if nbr_case % 10000 == 0:
      print(" calculus_integrate: generating case %d." % nbr_case)
    nbr_case += 1 

def generate_calculus_integrate_sample(vlist, ops, min_depth, max_depth,
                                       functions):
  """Randomly generate a symbolic integral dataset sample.

  Given an input expression, produce the indefinite integral.

  Args:
    vlist: Variable list. List of chars that can be used in the expression.
    ops: List of ExprOp instances. The allowed operators for the expression.
    min_depth: Expression trees will not have a smaller depth than this. 0 means
        there is just a variable. 1 means there is one operation.
    max_depth: Expression trees will not have a larger depth than this. To make
        all trees have the same depth, set this equal to `min_depth`.
    functions: Defines special functions. A dict mapping human readable string
        names, like "log", "exp", "sin", "cos", etc., to single chars. Each
        function gets a unique token, like "L" for "log".

  Returns:
    sample: String representation of the input. Will be of the form
        'var:expression'.
    target: String representation of the solution.
  """
  var_index = random.randrange(len(vlist))
  var = vlist[var_index]
  consts = vlist[:var_index] + vlist[var_index + 1:]

  depth = random.randrange(min_depth, max_depth + 1)
  expr = random_expr_with_required_var(depth, var, consts, ops)

  expr_str = str(expr)
  sample = var + ":" + expr_str
  target = format_sympy_expr(
      sympy.integrate(expr_str, sympy.Symbol(var)), functions=functions)
  return sample, target


# AlgebraConfig holds objects required to generate the algebra inverse
# dataset.
# vlist: Variable list. A list of chars.
# dlist: Numberical digit list. A list of chars.
# flist: List of special function names. A list of chars.
# functions: Dict of special function names. Maps human readable string names to
#     single char names used in flist.
# ops: Dict mapping op symbols (chars) to ExprOp instances.
# solve_ops: Encodes rules for how to algebraically cancel out each operation.
#     See doc-string for `algebra_inverse_solve`.
# int_encoder: Function that maps a string to a list of tokens. Use this to
#     encode an expression to feed into a model.
# int_decoder: Function that maps a list of tokens to a string. Use this to
#     convert model input or output into a human readable string.