Python sympy.lambdify() Examples

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(3, 1)

        x = sp.Matrix(sp.symbols('x y theta', real=True))
        u = sp.Matrix(sp.symbols('v w', real=True))

        f[0, 0] = u[0, 0] * sp.cos(x[2, 0])
        f[1, 0] = u[0, 0] * sp.sin(x[2, 0])
        f[2, 0] = u[1, 0]

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

def testSymbolicDims(self):
    p = builder.Base.Params()
    b = p.Instantiate()

    f1 = tshape.Shape(['kh', 'kw', 'idims', 'odims'])
    kh, kw, idims, odims = f1
    f2 = tshape.Shape([kh, kw, odims, odims])
    p = b._Seq('test', b._Conv2D('conv', f1, (2, 2)),
               b._Conv2D('conv', f2, (2, 2)), b._Bias('bias', odims))

    inp = tshape.Shape(['b', 'h', 'w', idims])
    b, h, w, _ = inp
    meta = p.cls.FPropMeta(p, inp)
    print('flops = ', meta.flops)
    out = meta.out_shapes[0]
    print('outputs = ', out)

    # sympy.lambdify can help us to do faster numerical evaluation.
    # Might be useful to build a "cost" model given a builder layer.
    f = sympy.lambdify([b, h, w, kh, kw, idims, odims], meta.flops, 'numpy')
    print('f.source = ', inspect.getsource(f))
    self.assertEqual(f(8, 224, 224, 3, 3, 8, 32), 925646848)
    self.assertEqual(f(8, 224, 224, 5, 5, 8, 32), 2569814016) 

def _is_non_decreasing(fn, q, bounds):
  """Verifies whether the function is non-decreasing within a range.

  Args:
    fn: Symbolic function of a single variable.
    q: The name of f's variable.
    bounds: Pair of (lower_bound, upper_bound) reals.

  Returns:
    True iff the function is non-decreasing in the range.
  """
  diff_fn = sp.diff(fn, q)  # Symbolically compute the derivative.
  diff_fn_lambdified = sp.lambdify(
      q,
      diff_fn,
      modules=[
          "numpy", {
              "erfc": scipy.special.erfc,
              "erfcinv": scipy.special.erfcinv
          }
      ])
  r = scipy.optimize.minimize_scalar(
      diff_fn_lambdified, bounds=bounds, method="bounded")
  assert r.success, "Minimizer failed to converge."
  return r.fun >= 0  # Check whether the derivative is non-negative. 

def _reset_samesymbols(self) -> None:
        """
        Set to a new function, assuming the same variables.
        """
        self.latex_repr = latex(self._symbolic_func)
        self._lambda_func = lambdify(
            self.symbols, self._symbolic_func) 

def convert_to_function(string: str, scale_by_k=False):
    """Using the sympy module, parse string input
    into a mathematical expression.
    Returns the original string, the latexified string,
    the mathematical expression in terms of sympy symbols,
    and a lambdified function
    """
    string = string.replace("^", "**")
    symbolic_function = parse_expr(string)
    if scale_by_k:
        latexstring = latex(symbolic_function*abc.k)
    else:
        latexstring = latex(symbolic_function)
    lambda_function = lambdify(abc.x, symbolic_function,
                               modules=module_list)
    string = string.replace('*', '')
    latexstring = "$" + latexstring + "$"
    return string, latexstring, \
           symbolic_function, lambda_function 

def expression_lambda(self) -> Callable:
        if self._expression_lambda is None:
            expression_lambda = sympy.lambdify(self.variables, self.underlying_expression,
                                                     [{'ceiling': ceiling}, 'numpy'])

            @functools.wraps(expression_lambda)
            def expression_wrapper(*args, **kwargs):
                result = expression_lambda(*args, **kwargs)
                if isinstance(result, sympy.NDimArray):
                    return numpy.array(result.tolist())
                elif isinstance(result, list):
                    return numpy.array(result).reshape(self.underlying_expression.shape)
                else:
                    return result.reshape(self.underlying_expression.shape)
            self._expression_lambda = expression_wrapper
        return self._expression_lambda 

def _is_non_decreasing(fn, q, bounds):
  """Verifies whether the function is non-decreasing within a range.

  Args:
    fn: Symbolic function of a single variable.
    q: The name of f's variable.
    bounds: Pair of (lower_bound, upper_bound) reals.

  Returns:
    True iff the function is non-decreasing in the range.
  """
  diff_fn = sp.diff(fn, q)  # Symbolically compute the derivative.
  diff_fn_lambdified = sp.lambdify(
      q,
      diff_fn,
      modules=[
          "numpy", {
              "erfc": scipy.special.erfc,
              "erfcinv": scipy.special.erfcinv
          }
      ])
  r = scipy.optimize.minimize_scalar(
      diff_fn_lambdified, bounds=bounds, method="bounded")
  assert r.success, "Minimizer failed to converge."
  return r.fun >= 0  # Check whether the derivative is non-negative. 

def set_symbolic_to_numeric_method(method): 
  ''' 
  Sets the method that all RBF instances will use for converting sympy
  expressions to numeric functions. This can be either "ufuncify" or
  "lambdify". "ufuncify" will write and compile C code for a numpy universal
  function, and "lambdify" will evaluate the sympy expression using
  python-level numpy functions. Calling this function will cause all caches of
  numeric functions to be cleared.
  '''
  global _SYMBOLIC_TO_NUMERIC_METHOD
  if method not in {'lambdify', 'ufuncify'}:
    raise ValueError(
      '`method` must be either "lambdify" or "ufuncify"')
            
  _SYMBOLIC_TO_NUMERIC_METHOD = method
  clear_rbf_caches()


## Instantiate some common RBFs
##################################################################### 

def vectors(self, mesh):
        h, Sig = self.fcts()

        f_hr = sympy.lambdify((r, t, z), h[0], 'numpy')
        f_ht = sympy.lambdify((r, t, z), h[1], 'numpy')
        f_hz = sympy.lambdify((r, t, z), h[2], 'numpy')

        f_sig = sympy.lambdify((r, t, z), Sig[0], 'numpy')

        hr = f_hr(mesh.gridEx[:, 0], mesh.gridEx[:, 1], mesh.gridEx[:, 2])
        ht = f_ht(mesh.gridEy[:, 0], mesh.gridEy[:, 1], mesh.gridEy[:, 2])
        hz = f_hz(mesh.gridEz[:, 0], mesh.gridEz[:, 1], mesh.gridEz[:, 2])

        sig = f_sig(mesh.gridCC[:, 0], mesh.gridCC[:, 1], mesh.gridCC[:, 2])

        return sig, np.r_[hr, ht, hz] 

def vectors(self, mesh):
        j, Sig = self.fcts()

        f_jr = sympy.lambdify((r, t, z), j[0], 'numpy')
        f_jt = sympy.lambdify((r, t, z), j[1], 'numpy')
        f_jz = sympy.lambdify((r, t, z), j[2], 'numpy')

        f_sig = sympy.lambdify((r, t, z), Sig[0], 'numpy')

        jr = f_jr(mesh.gridFx[:, 0], mesh.gridFx[:, 1], mesh.gridFx[:, 2])
        jt = f_jt(mesh.gridFy[:, 0], mesh.gridFy[:, 1], mesh.gridFy[:, 2])
        jz = f_jz(mesh.gridFz[:, 0], mesh.gridFz[:, 1], mesh.gridFz[:, 2])

        sig = f_sig(mesh.gridCC[:, 0], mesh.gridCC[:, 1], mesh.gridCC[:, 2])

        return sig, np.r_[jr, jt, jz] 

def _is_non_decreasing(fn, q, bounds):
  """Verifies whether the function is non-decreasing within a range.

  Args:
    fn: Symbolic function of a single variable.
    q: The name of f's variable.
    bounds: Pair of (lower_bound, upper_bound) reals.

  Returns:
    True iff the function is non-decreasing in the range.
  """
  diff_fn = sp.diff(fn, q)  # Symbolically compute the derivative.
  diff_fn_lambdified = sp.lambdify(
      q,
      diff_fn,
      modules=[
          "numpy", {
              "erfc": scipy.special.erfc,
              "erfcinv": scipy.special.erfcinv
          }
      ])
  r = scipy.optimize.minimize_scalar(
      diff_fn_lambdified, bounds=bounds, method="bounded")
  assert r.success, "Minimizer failed to converge."
  return r.fun >= 0  # Check whether the derivative is non-negative. 

def test_call_sim_with_args():
    a, a_unc, b, b_unc = 3.0, 0.1, 4.0, 0.1
    c = f_hypotenuse(a, b)
    m1 = PythagorasModel()
    data = {'PythagorasData': {'a': a, 'b': b, 'a_unc': a_unc, 'b_unc': b_unc}}
    m1.command('run', data=data)
    assert m1.registries['outputs']['c'].m == c
    assert m1.registries['outputs']['c'].u == UREG.cm
    x, y = sympy.symbols('x, y')
    z = sympy.sqrt(x * x + y * y)
    fx = sympy.lambdify((x, y), z.diff(x))
    fy = sympy.lambdify((x, y), z.diff(y))
    dz = np.sqrt(fx(a, b) ** 2 * a_unc ** 2 + fy(a, b) ** 2 * b_unc ** 2)
    c_unc = c * np.sqrt(m1.registries['outputs'].variance['c']['c'])
    LOGGER.debug('uncertainty in c is %g', c_unc)
    assert np.isclose(dz, c_unc.item())
    c_unc = c * m1.registries['outputs'].uncertainty['c']['c'].to('fraction')
    assert np.isclose(dz, c_unc.m.item())
    return m1 

def test_more_than_255_args():
    # SymPy's lambdify can handle at most 255 arguments
    # this is a proof of concept that this limitation does
    # not affect SymEngine's Lambdify class
    n = 257
    x = se.symarray('x', n)
    p, q, r = 17, 42, 13
    terms = [i*s for i, s in enumerate(x, p)]
    exprs = [se.add(*terms), r + x[0], -99]
    callback = se.Lambdify(x, exprs)
    input_arr = np.arange(q, q + n*n).reshape((n, n))
    out = callback(input_arr)
    ref = np.empty((n, 3))
    coeffs = np.arange(p, p + n, dtype=np.int64)
    for i in range(n):
        ref[i, 0] = coeffs.dot(np.arange(q + n*i, q + n*(i+1), dtype=np.int64))
        ref[i, 1] = q + n*i + r
    ref[:, 2] = -99
    assert np.allclose(out, ref) 

def _is_non_decreasing(fn, q, bounds):
  """Verifies whether the function is non-decreasing within a range.

  Args:
    fn: Symbolic function of a single variable.
    q: The name of f's variable.
    bounds: Pair of (lower_bound, upper_bound) reals.

  Returns:
    True iff the function is non-decreasing in the range.
  """
  diff_fn = sp.diff(fn, q)  # Symbolically compute the derivative.
  diff_fn_lambdified = sp.lambdify(
      q,
      diff_fn,
      modules=[
          "numpy", {
              "erfc": scipy.special.erfc,
              "erfcinv": scipy.special.erfcinv
          }
      ])
  r = scipy.optimize.minimize_scalar(
      diff_fn_lambdified, bounds=bounds, method="bounded")
  assert r.success, "Minimizer failed to converge."
  return r.fun >= 0  # Check whether the derivative is non-negative. 

def EvalExpr(value_type, x):
  """Evaluates x with symbol_to_value_map within the current context.

  Args:
    value_type: the target value type (see VALUE_TYPE).
    x: a sympy.Expr, an object, or a list/tuple of Exprs and objects.

  Returns:
    Evaluation result of 'x'.
  """
  if isinstance(x, (list, tuple)):
    return type(x)(EvalExpr(value_type, y) for y in x)
  elif isinstance(x, sympy.Expr):
    symbol_to_value_map = SymbolToValueMap.Get(value_type)
    if not symbol_to_value_map:
      return x
    # In theory the below should be equivalent to:
    #   y = x.subs(symbol_to_value_map).
    # In practice subs() doesn't work for when values are Tensors.
    k, v = list(zip(*(list(symbol_to_value_map.items()))))
    y = sympy.lambdify(k, x)(*v)
    return y
  else:
    return x 

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(6, 1)

        x = sp.Matrix(sp.symbols('rx ry vx vy t w', real=True))
        u = sp.Matrix(sp.symbols('gimbal T', real=True))

        f[0, 0] = x[2, 0]
        f[1, 0] = x[3, 0]
        f[2, 0] = 1 / self.m * sp.sin(x[4, 0] + u[0, 0]) * u[1, 0]
        f[3, 0] = 1 / self.m * (sp.cos(x[4, 0] + u[0, 0]) * u[1, 0] - self.m * self.g)
        f[4, 0] = x[5, 0]
        f[5, 0] = 1 / self.I * (-sp.sin(u[0, 0]) * u[1, 0] * self.r_T)

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

def evaluate_lambdified(expression: Union[sympy.Expr, numpy.ndarray],
                        variables: Sequence[str],
                        parameters: Dict[str, Union[numpy.ndarray, Number]],
                        lambdified) -> Tuple[Any, Any]:
    lambdified = lambdified or sympy.lambdify(variables, expression, _lambdify_modules)

    return lambdified(**parameters), lambdified 

def vectors(self, mesh):
        """ Get Vectors sig, sr. jx from sympy"""
        j, Sig = self.fcts()

        f_jr = sympy.lambdify((r, z), j[0], 'numpy')
        f_jz = sympy.lambdify((r, z), j[1], 'numpy')
        f_sigr = sympy.lambdify((r, z), Sig[0], 'numpy')
        f_sigz = sympy.lambdify((r, z), Sig[3], 'numpy')

        jr = f_jr(mesh.gridFx[:, 0], mesh.gridFx[:, 2])
        jz = f_jz(mesh.gridFz[:, 0], mesh.gridFz[:, 2])
        sigr = f_sigr(mesh.gridCC[:, 0], mesh.gridCC[:, 2])
        sigz = f_sigz(mesh.gridCC[:, 0], mesh.gridCC[:, 2])

        return np.c_[sigr, sigr, sigz], np.r_[jr, jz] 

def vectors(self, mesh):
        """ Get Vectors sig, sr. jx from sympy"""
        h, Sig = self.fcts()

        f_h = sympy.lambdify((r, z), h[0], 'numpy')
        f_sig = sympy.lambdify((r, z), Sig[0], 'numpy')

        ht = f_h(mesh.gridEy[:, 0], mesh.gridEy[:, 2])
        sig = f_sig(mesh.gridCC[:, 0], mesh.gridCC[:, 2])

        return sig, np.r_[ht] 

def evaluate_expression(
        expression: Union[sp.Expr, sp.Symbol], data: Mapping, data_override: Optional[Mapping] = None,
        function_mapping: Optional[Mapping[str, Callable]] = None) -> Array:
    """Evaluate a SymPy expression at data mapping variable names to arrays. Optionally, supplement the default suite of
    NumPy functions with a mapping from non-default function names to functions.
    """
    if expression.is_number:
        return np.asarray(float(expression), options.dtype)
    if expression.is_symbol:
        return get_symbol_data(expression, data, data_override)
    symbols = list(expression.free_symbols)
    modules = [function_mapping or {}, 'numpy']
    columns = (get_symbol_data(s, data, data_override) for s in symbols)
    return sp.lambdify(symbols, expression, modules)(*columns) 

def _numeric_solow_residual(self):
        """
        Vectorized, numpy-aware function defining the Solow residual.

        :getter: Return vectorized symbolic Solow residual.
        :type: function

        """
        if self.__numeric_solow_residual is None:
            tmp_args = [Y, K, L] + sym.symbols(list(self.params.keys()))
            self.__numeric_solow_residual = sym.lambdify(tmp_args,
                                                         self.solow_residual,
                                                         self._modules)
        return self.__numeric_solow_residual 

def lambdify_with_vector_args(args, expr, modules=DEFAULT_LAMBDIFY_MODULES):
    """
    A wrapper around sympy's lambdify where process_vector_args is used so
    generated callable can take arguments as either vector or individual
    components

    Parameters
    ----------
    args : list-like of sympy symbols
        Input arguments to the expression to call
    expr : sympy expression
        Expression to turn into a callable for numeric evaluation
    modules : list
        See lambdify documentation; passed directly as modules keyword.

    """
    new_args = process_vector_args(args)

    if sp.__version__ < '1.1' and hasattr(expr, '__len__'):
        expr = sp.Matrix(expr)

    f = sp.lambdify(new_args, expr, modules=modules)

    def lambda_function_with_vector_args(*func_args):
        new_func_args = process_vector_args(func_args)
        return np.array(f(*new_func_args))
    lambda_function_with_vector_args.__doc__ = f.__doc__
    return lambda_function_with_vector_args 

def _numeric_jacobian(self):
        """
        Vectorized, numpy-aware function defining the Jacobian matrix of
        partial derivatives.

        :getter: Return vectorized Jacobian matrix of partial derivatives.
        :type: function

        """
        if self.__numeric_jacobian is None:
            self.__numeric_jacobian = sym.lambdify(self._symbolic_args,
                                                   self._symbolic_jacobian,
                                                   self._modules)
        return self.__numeric_jacobian 

def execute_lambdify(ui, spacing=0.01, a=0.5, timesteps=500):
    """Execute diffusion stencil using vectorised numpy array accesses."""
    nx, ny = ui.shape
    dx2, dy2 = spacing**2, spacing**2
    dt = dx2 * dy2 / (2 * a * (dx2 + dy2))
    u = np.concatenate((ui, np.zeros_like(ui))).reshape((2, nx, ny))

    def diffusion_stencil():
        """Create stencil and substitutions for the diffusion equation"""
        p = sympy.Function('p')
        x, y, t, h, s = sympy.symbols('x y t h s')
        dx2 = p(x, y, t).diff(x, x).as_finite_difference([x - h, x, x + h])
        dy2 = p(x, y, t).diff(y, y).as_finite_difference([y - h, y, y + h])
        dt = p(x, y, t).diff(t).as_finite_difference([t, t + s])
        eqn = Eq(dt, a * (dx2 + dy2))
        stencil = solve(eqn, p(x, y, t + s))
        return stencil, (p(x, y, t), p(x + h, y, t), p(x - h, y, t),
                         p(x, y + h, t), p(x, y - h, t), s, h)
    stencil, subs = diffusion_stencil()
    kernel = sympy.lambdify(subs, stencil, 'numpy')

    # Execute timestepping loop with alternating buffers
    tstart = time.time()
    for ti in range(timesteps):
        t0 = ti % 2
        t1 = (ti + 1) % 2
        u[t1, 1:-1, 1:-1] = kernel(u[t0, 1:-1, 1:-1], u[t0, 2:, 1:-1],
                                   u[t0, :-2, 1:-1], u[t0, 1:-1, 2:],
                                   u[t0, 1:-1, :-2], dt, spacing)
    runtime = time.time() - tstart
    log("Lambdify: Diffusion with dx=%0.4f, dy=%0.4f, executed %d timesteps in %f seconds"
        % (spacing, spacing, timesteps, runtime))
    return u[ti % 2, :, :], runtime 

def vectors(self, mesh):
        """ Get Vectors sig, sr. jx from sympy"""
        j, Sig = self.fcts()

        f_jr = sympy.lambdify((r, z), j[0], 'numpy')
        f_jz = sympy.lambdify((r, z), j[1], 'numpy')
        f_sigr = sympy.lambdify((r, z), Sig[0], 'numpy')
        # f_sigz = sympy.lambdify((r,z), Sig[1], 'numpy')

        jr = f_jr(mesh.gridFx[:, 0], mesh.gridFx[:, 2])
        jz = f_jz(mesh.gridFz[:, 0], mesh.gridFz[:, 2])
        sigr = f_sigr(mesh.gridCC[:, 0], mesh.gridCC[:, 2])

        return sigr, np.r_[jr, jz] 

def _mpk(self):
        """
        :getter: Return vectorized symbolic marginal product capital.
        :type: function

        """
        if self.__mpk is None:
            args = [k] + sym.symbols(list(self.params.keys()))
            self.__mpk = sym.lambdify(args, self.marginal_product_capital,
                                      self._modules)
        return self.__mpk 

def delta(x: np.ndarray) -> np.ndarray:
    """
    Discrete approximation of the dirac delta function.
    """
    try:
        dx = (x[-1] - x[0])/len(x)
        return np.array([1e10 if (xi < (0. + dx/2) and xi > (0. - dx/2))
                         else 0. for xi in x])
    except:
        return 1e10 if (x < 0.01 and x > -0.01) \
               else 0.


# Dictionary of modules and user defined functions.
# Used for lambdify from sympy to parse input. 

def _intensive_output(self):
        """
        :getter: Return vectorized symbolic intensive aggregate production.
        :type: function

        """
        if self.__intensive_output is None:
            args = [k] + sym.symbols(list(self.params.keys()))
            self.__intensive_output = sym.lambdify(args, self.intensive_output,
                                                   self._modules)
        return self.__intensive_output 

def __init__(self, function_name: str,
                 param: Union[basic.Basic, str]) -> None:
        """
        The initializer. The parameter must be a
        string representation of a function, and it needs to
        be at least a function of x.
        """
        # Dictionary of modules and user defined functions.
        # Used for lambdify from sympy to parse input.
        if isinstance(param, str):
            param = parse_expr(param)
        if function_name == "x":
            function_name = "1.0*x"
        self._symbolic_func = parse_expr(function_name)
        symbol_set = self._symbolic_func.free_symbols
        if abc.k in symbol_set:
            k_param = parse_expr("k_param")
            self._symbolic_func = self._symbolic_func.subs(abc.k, k_param)
            symbol_set = self._symbolic_func.free_symbols
        symbol_list = list(symbol_set)
        if param not in symbol_list:
            raise VariableNotFoundError
        self.latex_repr = latex(self._symbolic_func)
        symbol_list.remove(param)
        self.parameters = symbol_list
        var_list = [param]
        var_list.extend(symbol_list)
        self.symbols = var_list
        self._lambda_func = lambdify(
            self.symbols, self._symbolic_func, modules=self.module_list) 

def global_iterator_to_indices(self, git=None):
        """
        Return sympy expressions translating global_iterator to loop indices.

        If global_iterator is given, an integer is returned
        """
        # unwind global iteration count into loop counters:
        base_loop_counters = {}
        global_iterator = symbol_pos_int('global_iterator')
        idiv = implemented_function(sympy.Function(str('idiv')), lambda x, y: x//y)
        total_length = 1
        last_incr = 1
        for var_name, start, end, incr in reversed(self._loop_stack):
            loop_var = symbol_pos_int(var_name)

            # This unspools the iterations:
            length = end-start  # FIXME is incr handled correct here?
            counter = start+(idiv(global_iterator*last_incr, total_length)*incr) % length
            total_length = total_length*length
            last_incr = incr

            base_loop_counters[loop_var] = sympy.lambdify(
                global_iterator,
                self.subs_consts(counter), modules=[numpy, {'Mod': numpy.mod}])

            if git is not None:
                try:  # Try to resolve to integer if global_iterator was given
                    base_loop_counters[loop_var] = sympy.Integer(self.subs_consts(counter))
                    continue
                except (ValueError, TypeError):
                    base_loop_counters[loop_var] = base_loop_counters[loop_var](git)

        return base_loop_counters