Python sympy.Function() Examples

def evaluate(self):
        # Average values if at a location not on the Function's grid
        if self._is_on_grid:
            return self
        weight = 1.0
        avg_list = [self]
        is_averaged = False
        for i, ir, d in zip(self.indices, self.indices_ref, self.dimensions):
            off = (i - ir)/d.spacing
            if not isinstance(off, sympy.Number) or int(off) == off:
                pass
            else:
                weight *= 1/2
                is_averaged = True
                avg_list = [(a.xreplace({i: i - d.spacing/2}) +
                             a.xreplace({i: i + d.spacing/2})) for a in avg_list]

        if not is_averaged:
            return self
        return weight * sum(avg_list) 

def testDiv(self):
    div = ops.Div(2, 3)
    self.assertEqual(str(div), '2/3')
    self.assertEqual(div.sympy(), sympy.Rational(2, 3))

    div = ops.Div(2, sympy.Rational(4, 5))
    self.assertEqual(str(div), '2/(4/5)')
    self.assertEqual(div.sympy(), sympy.Rational(5, 2))

    div = ops.Div(1, ops.Div(2, 3))
    self.assertEqual(str(div), '1/(2/3)')
    self.assertEqual(div.sympy(), sympy.Rational(3, 2))

    div = ops.Div(ops.Div(2, 3), 4)
    self.assertEqual(str(div), '(2/3)/4')
    self.assertEqual(div.sympy(), sympy.Rational(1, 6))

    div = ops.Div(2, ops.Mul(3, 4))
    self.assertEqual(str(div), '2/(3*4)')

    div = ops.Div(2, sympy.Function('f')(sympy.Symbol('x')))
    self.assertEqual(str(div), '2/f(x)') 

def __init__(self, *function_entities):
    """Initialize a `FunctionHandle`.

    Args:
      *function_entities: List of function letters and `Entity`s representing
          functions, to be composed.
    """
    self._functions = []
    for fn in function_entities:
      if isinstance(fn, str):
        functions = [sympy.Function(fn)]
      else:
        assert isinstance(fn, Entity)
        assert isinstance(fn.handle, FunctionHandle)
        functions = fn.handle.functions
      self._functions += functions 

def __init_finalize__(self, *args, **kwargs):
        super(SparseFunction, self).__init_finalize__(*args, **kwargs)
        self.interpolator = LinearInterpolator(self)
        # Set up sparse point coordinates
        coordinates = kwargs.get('coordinates', kwargs.get('coordinates_data'))
        if isinstance(coordinates, Function):
            self._coordinates = coordinates
        else:
            dimensions = (self.indices[-1], Dimension(name='d'))
            # Only retain the local data region
            if coordinates is not None:
                coordinates = np.array(coordinates)
            self._coordinates = SubFunction(name='%s_coords' % self.name, parent=self,
                                            dtype=self.dtype, dimensions=dimensions,
                                            shape=(self.npoint, self.grid.dim),
                                            space_order=0, initializer=coordinates,
                                            distributor=self._distributor)
            if self.npoint == 0:
                # This is a corner case -- we might get here, for example, when
                # running with MPI and some processes get 0-size arrays after
                # domain decomposition. We "touch" the data anyway to avoid the
                # case ``self._data is None``
                self.coordinates.data 

def as_tuple(item, type=None, length=None):
    """
    Force item to a tuple.

    Partly extracted from: https://github.com/OP2/PyOP2/.
    """
    # Empty list if we get passed None
    if item is None:
        t = ()
    elif isinstance(item, (str, sympy.Function)):
        t = (item,)
    else:
        # Convert iterable to list...
        try:
            t = tuple(item)
        # ... or create a list of a single item
        except (TypeError, NotImplementedError):
            t = (item,) * (length or 1)
    if length and not len(t) == length:
        raise ValueError("Tuple needs to be of length %d" % length)
    if type and not all(isinstance(i, type) for i in t):
        raise TypeError("Items need to be of type %s" % type)
    return t 

def _arg_check(self, args, intervals):
        """
        Check that ``args`` contains legal runtime values bound to ``self``.

        Raises
        ------
        InvalidArgument
            If, given the runtime values ``args``, an out-of-bounds array
            access would be performed, or if shape/dtype don't match with
            self's shape/dtype.
        """
        if self.name not in args:
            raise InvalidArgument("No runtime value for `%s`" % self.name)
        key = args[self.name]
        if len(key.shape) != self.ndim:
            raise InvalidArgument("Shape %s of runtime value `%s` does not match "
                                  "dimensions %s" %
                                  (key.shape, self.name, self.dimensions))
        if key.dtype != self.dtype:
            warning("Data type %s of runtime value `%s` does not match the "
                    "Function data type %s" % (key.dtype, self.name, self.dtype))
        for i, s in zip(self.dimensions, key.shape):
            i._arg_check(args, s, intervals[i]) 

def vars_as_functions(self):
        """
        :return: Turn the keys of this model into
            :class:`~sympy.core.function.Function`
            objects. This is done recursively so the chain rule can be applied
            correctly. This is done on the basis of `connectivity_mapping`.

            Example: for ``{y: a * x, z: y**2 + a}`` this returns
            ``{y: y(x, a), z: z(y(x, a), a)}``.
        """
        vars2functions = {}
        key = lambda arg: [isinstance(arg, Parameter), str(arg)]
        # Iterate over all symbols in this model in topological order, turning
        # each one into a function object recursively.
        for symbol in self.ordered_symbols:
            if symbol in self.connectivity_mapping:
                dependencies = self.connectivity_mapping[symbol]
                # Replace the dependency by it's function if possible
                dependencies = [vars2functions.get(dependency, dependency)
                               for dependency in dependencies]
                # sort by vars first, then params, and alphabetically within
                # each group
                dependencies = sorted(dependencies, key=key)
                vars2functions[symbol] = sympy.Function(symbol.name)(*dependencies)
        return vars2functions 

def local_indices(self):
        """
        Tuple of slices representing the global indices that logically
        belong to the calling MPI rank.

        Notes
        -----
        Given a Function ``f(x, y)`` with shape ``(nx, ny)``, when *not* using
        MPI this property will return ``(slice(0, nx-1), slice(0, ny-1))``. On
        the other hand, when MPI is used, the local ranges depend on the domain
        decomposition, which is carried by ``self.grid``.
        """
        if self._distributor is None:
            return tuple(slice(0, s) for s in self.shape)
        else:
            return tuple(self._distributor.glb_slices.get(d, slice(0, s))
                         for s, d in zip(self.shape, self.dimensions)) 

def test_conv11():
    x = sympy.Symbol("x")
    y = sympy.Symbol("y")
    x1 = Symbol("x")
    y1 = Symbol("y")
    f = sympy.Function("f")
    f1 = Function("f")

    e1 = diff(f(2*x, y), x)
    e2 = diff(f1(2*x1, y1), x1)
    e3 = diff(f1(2*x1, y1), y1)

    assert sympify(e1) == e2
    assert sympify(e1) != e3

    assert e2._sympy_() == e1
    assert e3._sympy_() != e1 

def test_conv10b():
    A = sympy.Matrix([[sympy.Symbol("x"), sympy.Symbol("y")],
                     [sympy.Symbol("z"), sympy.Symbol("t")]])
    assert sympify(A) == DenseMatrix(2, 2, [Symbol("x"), Symbol("y"),
                                            Symbol("z"), Symbol("t")])

    B = sympy.Matrix([[1, 2], [3, 4]])
    assert sympify(B) == DenseMatrix(2, 2, [Integer(1), Integer(2), Integer(3),
                                            Integer(4)])

    C = sympy.Matrix([[7, sympy.Symbol("y")],
                     [sympy.Function("g")(sympy.Symbol("z")), 3 + 2*sympy.I]])
    assert sympify(C) == DenseMatrix(2, 2, [Integer(7), Symbol("y"),
                                            function_symbol("g",
                                                            Symbol("z")),
                                            3 + 2*I]) 

def test_conv10():
    A = DenseMatrix(1, 4, [Integer(1), Integer(2), Integer(3), Integer(4)])
    assert (A._sympy_() == sympy.Matrix(1, 4,
                                        [sympy.Integer(1), sympy.Integer(2),
                                         sympy.Integer(3), sympy.Integer(4)]))

    B = DenseMatrix(4, 1, [Symbol("x"), Symbol("y"), Symbol("z"), Symbol("t")])
    assert (B._sympy_() == sympy.Matrix(4, 1,
                                        [sympy.Symbol("x"), sympy.Symbol("y"),
                                         sympy.Symbol("z"), sympy.Symbol("t")])
            )

    C = DenseMatrix(2, 2,
                    [Integer(5), Symbol("x"),
                     function_symbol("f", Symbol("x")), 1 + I])

    assert (C._sympy_() ==
            sympy.Matrix([[5, sympy.Symbol("x")],
                          [sympy.Function("f")(sympy.Symbol("x")),
                           1 + sympy.I]])) 

def _coeff_symbol(self):
        if self.coefficients == 'symbolic':
            return sympy.Function('W')
        else:
            raise ValueError("Function was not declared with symbolic "
                             "coefficients.") 

def __init__(self, **traits):
        TaylorPoly.__init__(self, **traits)
        #Declare the analytical function
  
        Z=sympy.Function("Z")
        Ax,Ay,Kx,Ky=sympy.symbols(("Ax","Ay","Kx","Ky"))
        x, y =sympy.symbols('xy')
        Z=(Ax*x**2+Ay*y**2)/(1+sympy.sqrt(1-(1+Kx)*Ax**2*x**2-(1+Ky)*Ay**2*y**2));
        
        #Calculate taylor polynomial coheficients
        cohef=[[Z, ],]
        order=self.n
        for i in range(0, order+1, 2):
            if i!=0:
                cohef.append([sympy.diff(cohef[i/2-1][0], y, 2), ])
            for j in range(2, order-i+1, 2):
                cohef[i/2].append(sympy.diff(cohef[i/2][j/2 -1], x, 2))
        

        A_x=self.Ax
        A_y=self.Ay
        K_x=self.Kx
        K_y=self.Ky
        
        c=zeros((self.n+1, self.n+1))
        for i in range(0, order/2+1):
            for j in range(0,order/2- i+1):
                cohef[j][i]=cohef[j][i].subs(x, 0).subs(y, 0).subs(Ax, A_x).subs(Ay, A_y).subs(Kx, K_x).subs(Ky, K_y)/(sympy.factorial(2*i)*sympy.factorial(2*j))
                c[2*j, 2*i]=cohef[j][i].evalf()
        
        # Add the high order corrections
        if len(self.ho_cohef.shape)==2:
            cx, cy = c.shape 
            dx, dy =self.ho_cohef.shape
            mx=array((cx, dx)).max()
            my=array((cy, dy)).max()
            self.cohef=zeros((mx, my))
            self.cohef[0:cx, 0:cy]=c
            self.cohef[0:dy, 0:dy]=self.cohef[0:dy, 0:dy]+self.ho_cohef
        else:
            self.cohef=c 

def visit_NumpyAttr(self, node):
        if node.name in SYM_UNARY_FUNCTIONS:
            return SYM_UNARY_FUNCTIONS[node.name]
        else:
            return sp.Function('np.' + node.name)
    # TODO: implement this!
    # def visit_NumpyContraction(self, node): ??? 

def __init__(self):
        for name in list(SYM_UNARY_FUNCTIONS.keys()):
            if not name in NPY_UNARY_FUNCTIONS and not name in NPY_CAST_FUNCTIONS:
                raise Exception("Unknown Function {0}".format(name))
            setattr(self, "visit_{0}".format(name), self.npUnaryFunction_visit) 

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 __new__(cls, name, arguments=None):
        arguments = as_tuple(arguments)
        obj = Function.__new__(cls, name, *arguments)
        obj._name = name
        obj._arguments = arguments
        return obj 

def __new__(cls, *args, **kwargs):
        options = kwargs.get('options', {})

        key = cls._cache_key(*args, **kwargs)
        obj = cls._cache_get(key)

        if obj is not None:
            newobj = sympy.Matrix.__new__(cls, *args, **options)
            newobj.__init_cached__(key)
            return newobj

        name = kwargs.get('name')
        indices, _ = cls.__indices_setup__(**kwargs)

        # Create new, unique type instance from cls and the symbol name
        newcls = type(name, (cls,), dict(cls.__dict__))

        # Create the new Function object and invoke __init__
        comps = cls.__subfunc_setup__(*args, **kwargs)
        newobj = sympy.ImmutableDenseMatrix.__new__(newcls, comps)
        # Initialization. The following attributes must be available
        newobj._indices = indices
        newobj._name = name
        newobj._dtype = cls.__dtype_setup__(**kwargs)
        newobj.__init_finalize__(*args, **kwargs)

        # Store new instance in symbol cache
        Cached.__init__(newobj, newcls)

        return newobj 

def __indices_setup__(cls, **kwargs):
        dimensions = kwargs.get('dimensions')
        staggered = kwargs.get('staggered')
        if dimensions is None:
            save = kwargs.get('save')
            grid = kwargs.get('grid')
            time_dim = kwargs.get('time_dim')

            if time_dim is None:
                time_dim = grid.time_dim if isinstance(save, int) else grid.stepping_dim
            elif not (isinstance(time_dim, Dimension) and time_dim.is_Time):
                raise TypeError("`time_dim` must be a time dimension")
            dimensions = list(Function.__indices_setup__(**kwargs)[0])
            dimensions.insert(cls._time_position, time_dim)
        return Function.__indices_setup__(dimensions=dimensions, staggered=staggered) 

def _coordinate_indices(self):
        """Symbol for each grid index according to the coordinates."""
        indices = self.grid.dimensions
        return tuple([INT(sympy.Function('floor')((c - o) / i.spacing))
                      for c, o, i in zip(self._coordinate_symbols, self.grid.origin,
                                         indices[:self.grid.dim])]) 

def inject(self, field, expr, offset=0, u_t=None, p_t=None):
        """
        Generate equations injecting an arbitrary expression into a field.

        Parameters
        ----------
        field : Function
            Input field into which the injection is performed.
        expr : expr-like
            Injected expression.
        offset : int, optional
            Additional offset from the boundary.
        u_t : expr-like, optional
            Time index at which the interpolation is performed.
        p_t : expr-like, optional
            Time index at which the result of the interpolation is stored.
        """
        # Apply optional time symbol substitutions to field and expr
        if u_t is not None:
            field = field.subs({field.time_dim: u_t})
        if p_t is not None:
            expr = expr.subs({self.time_dim: p_t})

        return super(SparseTimeFunction, self).inject(field, expr, offset=offset)

    # Pickling support 

def eval(self, string: str, **_: Any) -> Array:
        """Parse a SymPy expression from a string and evaluate it at data represented as the environment's only
        namespace.
        """
        data = self._namespaces[0].copy()

        # parse the SymPy expression, preserving the function that marks variables as categorical
        expression = parse_expression(string, mark_categorical=True)

        # replace categorical variables with unicode objects and explicitly mark them as categorical so that labels are
        #   unique and so that all categorical variables are treated the same
        C = sp.Function('C')
        for symbol in expression.free_symbols:
            if not np.issubdtype(data[symbol.name].dtype, getattr(np, 'number')):
                expression = expression.replace(symbol, C(symbol))
                data[symbol.name] = data[symbol.name].astype(np.unicode_).astype(np.object_)

        # evaluate the expression and handle universally-marked categorical variables with a non-default coding class
        evaluated = evaluate_expression(expression, self._namespaces[0], function_mapping={
            'C': functools.partial(patsy.builtins.C, contrast=CategoricalTreatment)
        })

        # if the evaluated expression is a scalar, it is a constant that needs to be repeated
        if isinstance(evaluated, (numbers.Number, np.ndarray)) and np.asarray(evaluated).size == 1:
            size = next(iter(data.values())).shape[0]
            evaluated = np.ones(size) * evaluated
        return evaluated 

def test_conv8b():
    e1 = sympy.Function("f")(sympy.Symbol("x"))
    e2 = sympy.Function("g")(sympy.Symbol("x"), sympy.Symbol("y"))
    assert sympify(e1) == function_symbol("f", Symbol("x"))
    assert sympify(e2) != function_symbol("f", Symbol("x"))
    assert sympify(e2) == function_symbol("g", Symbol("x"), Symbol("y"))

    e3 = sympy.Function("q")(sympy.Symbol("t"))
    assert sympify(e3) == function_symbol("q", Symbol("t"))
    assert sympify(e3) != function_symbol("f", Symbol("t"))
    assert sympify(e3) != function_symbol("q", Symbol("t"), Symbol("t")) 

def test_conv8():
    e1 = function_symbol("f", Symbol("x"))
    e2 = function_symbol("g", Symbol("x"), Symbol("y"))
    assert e1._sympy_() == sympy.Function("f")(sympy.Symbol("x"))
    assert e2._sympy_() != sympy.Function("f")(sympy.Symbol("x"))
    assert (e2._sympy_() ==
            sympy.Function("g")(sympy.Symbol("x"), sympy.Symbol("y")))

    e3 = function_symbol("q", Symbol("t"))
    assert e3._sympy_() == sympy.Function("q")(sympy.Symbol("t"))
    assert e3._sympy_() != sympy.Function("f")(sympy.Symbol("t"))
    assert (e3._sympy_() !=
            sympy.Function("q")(sympy.Symbol("t"), sympy.Symbol("t"))) 

def __init__(self, base=None, fct=None, deriv=None, on=True, debug=False):
        if on:
            OS = 'unix'
            if 'win' in sys.platform and 'darwin' not in sys.platform:
                OS = 'win'

            if base is None:
                Eprint.base = Eprint.ColorCode[Eprint.defaults[(OS, 'base')]]
            else:
                Eprint.base = Eprint.ColorCode[base]
            if fct is None:
                Eprint.fct = Eprint.ColorCode[Eprint.defaults[(OS, 'fct')]]
            else:
                Eprint.fct = Eprint.ColorCode[fct]
            if deriv is None:
                Eprint.deriv = Eprint.ColorCode[Eprint.defaults[(OS, 'deriv')]]
            else:
                Eprint.deriv = Eprint.ColorCode[deriv]
            Eprint.normal = '\033[0m'

            if debug:
                print('Enhanced Printing is on:')
                print('Base/Blade color is ' + Eprint.InvColorCode[Eprint.base])
                print('Function color is ' + Eprint.InvColorCode[Eprint.fct])
                print('Derivative color is ' + Eprint.InvColorCode[Eprint.deriv] + '\n')

            Eprint.base = '\033[' + Eprint.base + 'm'
            Eprint.fct = '\033[' + Eprint.fct + 'm'
            Eprint.deriv = '\033[' + Eprint.deriv + 'm' 

def find_functions(expr):
    f_lst = []
    for f in list(expr.atoms(Function)):
        if str(f) not in GaPrinter.function_names:
            f_lst.append(f)
    f_lst += list(expr.atoms(Derivative))
    return f_lst 

def test_not_fortran():
    x = symbols('x')
    g = Function('g')
    assert fcode(
        gamma(x)) == "C     Not Fortran:\nC     gamma(x)\n      gamma(x)"
    assert fcode(Integral(sin(x))) == "C     Not Fortran:\nC     Integral(sin(x), x)\n      Integral(sin(x), x)"
    assert fcode(g(x)) == "C     Not Fortran:\nC     g(x)\n      g(x)" 

def test_printmethod():
    x = symbols('x')

    class nint(Function):
        def _fcode(self, printer):
            return "nint(%s)" % printer._print(self.args[0])
    assert fcode(nint(x)) == "      nint(x)" 

def python(expr, **settings):
    """Return Python interpretation of passed expression
    (can be passed to the exec() function without any modifications)"""

    printer = PythonPrinter(settings)
    exprp = printer.doprint(expr)

    result = ''
    # Returning found symbols and functions
    renamings = {}
    for symbolname in printer.symbols:
        newsymbolname = symbolname
        # Escape symbol names that are reserved python keywords
        if kw.iskeyword(newsymbolname):
            while True:
                newsymbolname += "_"
                if (newsymbolname not in printer.symbols and
                        newsymbolname not in printer.functions):
                    renamings[sympy.Symbol(
                        symbolname)] = sympy.Symbol(newsymbolname)
                    break
        result += newsymbolname + ' = Symbol(\'' + symbolname + '\')\n'

    for functionname in printer.functions:
        newfunctionname = functionname
        # Escape function names that are reserved python keywords
        if kw.iskeyword(newfunctionname):
            while True:
                newfunctionname += "_"
                if (newfunctionname not in printer.symbols and
                        newfunctionname not in printer.functions):
                    renamings[sympy.Function(
                        functionname)] = sympy.Function(newfunctionname)
                    break
        result += newfunctionname + ' = Function(\'' + functionname + '\')\n'

    if not len(renamings) == 0:
        exprp = expr.subs(renamings)
    result += 'e = ' + printer._str(exprp)
    return result 

def _print_expint(self, e):
        from sympy import Function
        if e.args[0].is_Integer and self._use_unicode:
            return self._print_Function(Function('E_%s' % e.args[0])(e.args[1]))
        return self._print_Function(e)