Python sympy.Piecewise() Examples

def test_PiecewisePoly__sympy():
    import sympy as sp
    Poly = mk_Poly('T')
    p1 = Poly([0, 1, 0.1])
    p2 = Poly([0, 3, -.1])

    TPiecewisePoly = mk_PiecewisePoly('temperature')
    tpwp = TPiecewisePoly([2, 2, 0, 10, 2, 10, 20, 0, 1, 0.1, 0, 3, -.1])
    x = sp.Symbol('x')
    res = tpwp.eval_poly({'temperature': x}, backend=sp)
    assert isinstance(res, sp.Piecewise)
    assert res.args[0][0] == 1 + 0.1*x
    assert res.args[0][1] == sp.And(0 <= x, x <= 10)
    assert res.args[1][0] == 3 - 0.1*x
    assert res.args[1][1] == sp.And(10 <= x, x <= 20)

    with pytest.raises(ValueError):
        tpwp.from_polynomials([(0, 10), (10, 20)], [p1, p2]) 

def duration(self) -> ExpressionScalar:
        step_size = self._loop_range.step.sympified_expression
        loop_index = sympy.symbols(self._loop_index)
        sum_index = sympy.symbols(self._loop_index)

        # replace loop_index with sum_index dependable expression
        body_duration = self.body.duration.sympified_expression.subs({loop_index: self._loop_range.start.sympified_expression + sum_index*step_size})

        # number of sum contributions
        step_count = sympy.ceiling((self._loop_range.stop.sympified_expression-self._loop_range.start.sympified_expression) / step_size)
        sum_start = 0
        sum_stop = sum_start + (sympy.functions.Max(step_count, 1) - 1)

        # expression used if step_count >= 0
        finite_duration_expression = sympy.Sum(body_duration, (sum_index, sum_start, sum_stop))

        duration_expression = sympy.Piecewise((0, step_count <= 0),
                                              (finite_duration_expression, True))

        return ExpressionScalar(duration_expression) 

def _print_Piecewise(self, expr, **kwargs):
        import numpy as np
        e, cond = expr.args[0].args
        if len(expr.args) == 1:
            return tt.switch(self._print(cond, **kwargs),
                             self._print(e, **kwargs),
                             np.nan)
        return tt.switch(self._print(cond, **kwargs),
                         self._print(e, **kwargs),
                         self._print(sympy.Piecewise(*expr.args[1:]), **kwargs)) 

def _needs_mul_brackets(self, expr, last=False):
        """
        Returns True if the expression needs to be wrapped in brackets when
        printed as part of a Mul, False otherwise. This is True for Add,
        but also for some container objects that would not need brackets
        when appearing last in a Mul, e.g. an Integral. ``last=True``
        specifies that this expr is the last to appear in a Mul.
        """
        from sympy import Integral, Piecewise, Product, Sum
        return expr.is_Add or (not last and
            any([expr.has(x) for x in (Integral, Piecewise, Product, Sum)])) 

def test_logic():
    x = true
    y = false
    x1 = sympy.true
    y1 = sympy.false

    assert And(x, y) == And(x1, y1)
    assert And(x1, y) == And(x1, y1)
    assert And(x, y)._sympy_() == sympy.And(x1, y1)
    assert sympify(sympy.And(x1, y1)) == And(x, y)

    assert Or(x, y) == Or(x1, y1)
    assert Or(x1, y) == Or(x1, y1)
    assert Or(x, y)._sympy_() == sympy.Or(x1, y1)
    assert sympify(sympy.Or(x1, y1)) == Or(x, y)

    assert Not(x) == Not(x1)
    assert Not(x1) == Not(x1)
    assert Not(x)._sympy_() == sympy.Not(x1)
    assert sympify(sympy.Not(x1)) == Not(x)

    assert Xor(x, y) == Xor(x1, y1)
    assert Xor(x1, y) == Xor(x1, y1)
    assert Xor(x, y)._sympy_() == sympy.Xor(x1, y1)
    assert sympify(sympy.Xor(x1, y1)) == Xor(x, y)

    x = Symbol("x")
    x1 = sympy.Symbol("x")

    assert Piecewise((x, x < 1), (0, True)) == Piecewise((x1, x1 < 1), (0, True))
    assert Piecewise((x, x1 < 1), (0, True)) == Piecewise((x1, x1 < 1), (0, True))
    assert Piecewise((x, x < 1), (0, True))._sympy_() == sympy.Piecewise((x1, x1 < 1), (0, True))
    assert sympify(sympy.Piecewise((x1, x1 < 1), (0, True))) == Piecewise((x, x < 1), (0, True))

    assert Contains(x, Interval(1, 1)) == Contains(x1, Interval(1, 1))
    assert Contains(x, Interval(1, 1))._sympy_() == sympy.Contains(x1, Interval(1, 1))
    assert sympify(sympy.Contains(x1, Interval(1, 1))) == Contains(x, Interval(1, 1)) 

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 f(self, yvec, params):
        import sympy as sp
        f_penalty = [sp.Piecewise((yi**2, yi < 0), (0, True)) for yi in yvec]
        return super(_NumSysLinNegPenalty, self).f(yvec, params) + f_penalty 

def f(self, yvec, params):
        import sympy
        ymax = self.eqsys.upper_conc_bounds(
            params[:self.eqsys.ns],
            min_=lambda a, b: sympy.Piecewise((a, a < b), (b, True)))
        ytanh = [yimax*(4 + 5*sympy.tanh(yi))/8
                 for yimax, yi in zip(ymax, yvec)]
        return NumSysLin.f(self, ytanh, params) 

def test_PiecewiseTPolyMassAction__sympy():
    import sympy as sp
    tp1 = TPoly([10, 0.1])
    tp2 = ShiftedTPoly([273.15, 37.315, -0.1])
    pwp = MassAction(TPiecewise([0, tp1, 273.15, tp2, 373.15]))
    T = sp.Symbol('T')
    r = Reaction({'A': 2, 'B': 1}, {'C': 1}, inact_reac={'B': 1})
    res1 = pwp({'A': 11, 'B': 13, 'temperature': T}, backend=sp, reaction=r)
    ref1 = 11**2 * 13 * sp.Piecewise(
        (10+0.1*T, sp.And(0 <= T, T <= 273.15)),
        (37.315 - 0.1*(T-273.15), sp.And(273.15 <= T, T <= 373.15)),
        (sp.Symbol('NAN'), True)
    )
    assert res1 == ref1 

def test_create_Piecewise_Poly__sympy():
    import sympy as sp
    Poly = create_Poly('Tmpr')
    p1 = Poly([1, 0.1])
    p2 = Poly([3, -.1])

    TPw = create_Piecewise('Tmpr')
    pw = TPw([0, p1, 10, p2, 20])
    x = sp.Symbol('x')
    res = pw({'Tmpr': x}, backend=sp)
    assert isinstance(res, sp.Piecewise)
    assert res.args[0][0] == 1 + 0.1*x
    assert res.args[0][1] == sp.And(0 <= x, x <= 10)
    assert res.args[1][0] == 3 - 0.1*x
    assert res.args[1][1] == sp.And(10 <= x, x <= 20) 

def _print_hyper(self, e):
        # FIXME refactor Matrix, Piecewise, and this into a tabular environment
        ap = [self._print(a) for a in e.ap]
        bq = [self._print(b) for b in e.bq]

        P = self._print(e.argument)
        P.baseline = P.height()//2

        # Drawing result - first create the ap, bq vectors
        D = None
        for v in [ap, bq]:
            D_row = self._hprint_vec(v)
            if D is None:
                D = D_row       # first row in a picture
            else:
                D = prettyForm(*D.below(' '))
                D = prettyForm(*D.below(D_row))

        # make sure that the argument `z' is centred vertically
        D.baseline = D.height()//2

        # insert horizontal separator
        P = prettyForm(*P.left(' '))
        D = prettyForm(*D.right(' '))

        # insert separating `|`
        D = self._hprint_vseparator(D, P)

        # add parens
        D = prettyForm(*D.parens('(', ')'))

        # create the F symbol
        above = D.height()//2 - 1
        below = D.height() - above - 1

        sz, t, b, add, img = annotated('F')
        F = prettyForm('\n' * (above - t) + img + '\n' * (below - b),
                       baseline=above + sz)
        add = (sz + 1)//2

        F = prettyForm(*F.left(self._print(len(e.ap))))
        F = prettyForm(*F.right(self._print(len(e.bq))))
        F.baseline = above + add

        D = prettyForm(*F.right(' ', D))

        return D 

def _print_Mul(self, product):
        a = []  # items in the numerator
        b = []  # items that are in the denominator (if any)

        if self.order not in ('old', 'none'):
            args = product.as_ordered_factors()
        else:
            args = product.args

        # Gather terms for numerator/denominator
        for item in args:
            if item.is_commutative and item.is_Pow and item.exp.is_Rational and item.exp.is_negative:
                b.append(C.Pow(item.base, -item.exp))
            elif item.is_Rational and item is not S.Infinity:
                if item.p != 1:
                    a.append( C.Rational(item.p) )
                if item.q != 1:
                    b.append( C.Rational(item.q) )
            else:
                a.append(item)

        from sympy import Integral, Piecewise, Product, Sum

        # Convert to pretty forms. Add parens to Add instances if there
        # is more than one term in the numer/denom
        for i in xrange(0, len(a)):
            if (a[i].is_Add and len(a) > 1) or (i != len(a) - 1 and
                    isinstance(a[i], (Integral, Piecewise, Product, Sum))):
                a[i] = prettyForm(*self._print(a[i]).parens())
            else:
                a[i] = self._print(a[i])

        for i in xrange(0, len(b)):
            if (b[i].is_Add and len(b) > 1) or (i != len(b) - 1 and
                    isinstance(b[i], (Integral, Piecewise, Product, Sum))):
                b[i] = prettyForm(*self._print(b[i]).parens())
            else:
                b[i] = self._print(b[i])

        # Construct a pretty form
        if len(b) == 0:
            return prettyForm.__mul__(*a)
        else:
            if len(a) == 0:
                a.append( self._print(S.One) )
            return prettyForm.__mul__(*a)/prettyForm.__mul__(*b)

    # A helper function for _print_Pow to print x**(1/n) 

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 _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) 

def _add_diff_to_cache(self, diff):
    '''     
    Symbolically differentiates the RBF and then converts the expression to a
    function which can be evaluated numerically.
    '''   
    logger.debug('Creating a numerical function for the RBF %s with '
                 'the derivative %s ...' % (self, str(diff)))
    dim = len(diff)
    c_sym = sympy.symbols('c:%s' % dim)
    x_sym = sympy.symbols('x:%s' % dim)    
    r_sym = sympy.sqrt(sum((xi-ci)**2 for xi, ci in zip(x_sym, c_sym)))

    # substitute 'r' in the RBF expression with the cartesian spatial variables
    # and differentiate the RBF with respect to them
    expr = self.expr.subs(_R, r_sym)            
    for xi, order in zip(x_sym, diff):
      if order == 0:
        continue

      expr = expr.diff(*(xi,)*order)

    # if `tol` is given, form a separate expression for the RBF near its center
    if self.tol is not None:
      if diff in self.limits:
        # use a user-specified limit if available      
        lim = self.limits[diff]

      else: 
        logger.debug('Symbolically evaluating the RBF at its center ...')
        # evaluate the limit of the RBF at (x0=tol+c0, x1=c1, x2=c2, ...) as
        # tol goes to zero.
        lim = expr.subs(zip(x_sym[1:], c_sym[1:]))
        lim = lim.simplify()
        lim = lim.limit(x_sym[0], c_sym[0])
        logger.debug('Value of the RBF at its center: %s' % lim)

      # create a piecewise symbolic function which is `lim` when `r_sym < tol`
      # and `expr` otherwise
      expr = sympy.Piecewise((lim, r_sym < self.tol), (expr, True)) 

    if _SYMBOLIC_TO_NUMERIC_METHOD == 'ufuncify':      
      func = ufuncify(x_sym + c_sym + (_EPS,), expr, backend='numpy')

    elif _SYMBOLIC_TO_NUMERIC_METHOD == 'lambdify':
      func = lambdify(x_sym + c_sym + (_EPS,), expr, modules=['numpy'])

    else:
      raise ValueError()          
        
    self._cache[diff] = func
    logger.debug('The numeric function has been created and cached')