如何使用 numpy 在 Python 中定义分段函数？

Question

以下是我想在 python 中实现的功能。定义函数时出现类型错误。我尝试使用numpy.piecewise 函数对象进行定义，并且仅使用elif 命令作为定义。然后我希望能够在不同的点以及 f(X-1) 等表达式评估这个函数

这是我的代码：

from numpy import piecewise 
from scipy import *
from sympy.abc import x
from sympy.utilities.lambdify import lambdify, implemented_function
from sympy import Function
from sympy import *
h = 0.5 
a = -1
n = 2
x = Symbol('x')
expr = piecewise((0, x-a <=  -2*h), ((1/6)*(2*h+(x-a))**3, -2*h<=x-a<=-h), (2*h**3/3-0.5*(x-a)**2*(2*h+(x-a)), -h<= x-a<= 0), (2*(h**3/3)-0.5*(x-a)**2*(2*h+(x-a)), 0<=x-a<=2*h), ((1/6)*(2*h-(x-a))**3, h<=x-a<=2*h), (0, x-a<=2*h))
p = lambdify((x, a,b,h), expr)

def basis(x,a,b, h):
    if x <= a-2*h:
        return 0;
    elif (x<=a-h) or (x >=2*h):
        return (1/6)*(2*h+(x-a))**3
    elif  (x-a<= 0) or (x-a >= -h):
        return (2*h**3/3-0.5*(x-a)**2*(2*h+(x-a)));
    elif (x<=2*h+a) or (x >= 0):
        return  (2*(h**3/3)-0.5*(x-a)**2*(2*h+(x-a)));
    elif (x<=a+2*h) or (x >= h):
        return (1/6)*(2*h-(x-a))**3; 
    elif x-a<=2*h:
        return 0

basis(x, -1,0.5,0)

我得到这个的两种方式：

raise TypeError("cannot determine truth value of Relational")

TypeError: cannot determine truth value of Relational

Answer 1

您可以使用 sympy 的 lambdify 函数来生成 numpy 分段函数。这是一个更简单的示例，但显示了总体思路：

In [15]: from sympy import symbols, Piecewise                                                                                               

In [16]: x, a = symbols('x, a')                                                                                                   

In [17]: expr = Piecewise((x, x>a), (0, True))                                                                                    

In [18]: expr                                                                                                                     
Out[18]: 
⎧x  for a < x
⎨            
⎩0  otherwise

In [19]: from sympy import lambdify                                                                                               

In [20]: fun = lambdify((x, a), expr)                                                                                             

In [21]: fun([1, 3], [4, 2])                                                                                                      
Out[21]: array([0., 3.])

In [22]: import inspect                                                                                                           

In [23]: print(inspect.getsource(fun))                                                                                            
def _lambdifygenerated(x, a):
    return (select([less(a, x),True], [x,0], default=nan))

Answer 2

抱歉此答案的长度，但我认为您需要查看完整的调试过程。我不得不查看回溯并测试小块代码以确定确切的问题。我见过很多numpy ambiguity错误，但不是这个sympy关系错误。

===

让我们看一下整个回溯，而不仅仅是一行。至少我们需要识别您的代码的哪一行正在产生问题。

In [4]: expr = np.piecewise((0, x-a <=  -2*h), ((1/6)*(2*h+(x-a))**3, -2*h<=x-a<
   ...: =-h), (2*h**3/3-0.5*(x-a)**2*(2*h+(x-a)), -h<= x-a<= 0), (2*(h**3/3)-0.5
   ...: *(x-a)**2*(2*h+(x-a)), 0<=x-a<=2*h), ((1/6)*(2*h-(x-a))**3, h<=x-a<=2*h)
   ...: , (0, x-a<=2*h))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-4-893bb4b36321> in <module>
----> 1 expr = np.piecewise((0, x-a <=  -2*h), ((1/6)*(2*h+(x-a))**3, -2*h<=x-a<=-h), (2*h**3/3-0.5*(x-a)**2*(2*h+(x-a)), -h<= x-a<= 0), (2*(h**3/3)-0.5*(x-a)**2*(2*h+(x-a)), 0<=x-a<=2*h), ((1/6)*(2*h-(x-a))**3, h<=x-a<=2*h), (0, x-a<=2*h))

/usr/local/lib/python3.8/dist-packages/sympy/core/relational.py in __nonzero__(self)
    382 
    383     def __nonzero__(self):
--> 384         raise TypeError("cannot determine truth value of Relational")
    385 
    386     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

虽然np.piecewise 是一个numpy 函数，但因为x 是一个sympy.Symbol，所以方程是sympy 表达式。 numpy和sympy @ and and leation集成。有些东西有效，而其他许多则无效。

你试过一个小表情吗？良好的编程练习是从小块开始，确保这些工作首先。

让我们尝试更小的东西：

In [8]: expr = np.piecewise((0, x-a <=  -2*h),
   ...:  ((1/6)*(2*h+(x-a))**3, -2*h<=x-a<=-h))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-8-37ff62e49efb> in <module>
      1 expr = np.piecewise((0, x-a <=  -2*h),
----> 2  ((1/6)*(2*h+(x-a))**3, -2*h<=x-a<=-h))

/usr/local/lib/python3.8/dist-packages/sympy/core/relational.py in __nonzero__(self)
    382 
    383     def __nonzero__(self):
--> 384         raise TypeError("cannot determine truth value of Relational")
    385 
    386     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

和较小的碎片：

In [10]: (0, x-a <=  -2*h)
Out[10]: (0, x + 1 ≤ -1.0)

In [11]: ((1/6)*(2*h+(x-a))**3, -2*h<=x-a<=-h)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-7bd9f95d077d> in <module>
----> 1 ((1/6)*(2*h+(x-a))**3, -2*h<=x-a<=-h)

/usr/local/lib/python3.8/dist-packages/sympy/core/relational.py in __nonzero__(self)
    382 
    383     def __nonzero__(self):
--> 384         raise TypeError("cannot determine truth value of Relational")
    385 
    386     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

In [12]: (1/6)*(2*h+(x-a))**3
Out[12]: 
                            3
1.33333333333333⋅(0.5⋅x + 1) 

但是：

In [13]: -2*h<=x-a<=-h
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-13-5ffb419cd443> in <module>
----> 1 -2*h<=x-a<=-h

/usr/local/lib/python3.8/dist-packages/sympy/core/relational.py in __nonzero__(self)
    382 
    383     def __nonzero__(self):
--> 384         raise TypeError("cannot determine truth value of Relational")
    385 
    386     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

简化进一步：

In [14]: 0 < x < 3
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-59ba4ce00627> in <module>
----> 1 0 < x < 3

/usr/local/lib/python3.8/dist-packages/sympy/core/relational.py in __nonzero__(self)
    382 
    383     def __nonzero__(self):
--> 384         raise TypeError("cannot determine truth value of Relational")
    385 
    386     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

a < b < c @允许常规Python变量和标量，它不适用于numpy阵列，并且显然不适用于sympy变量。

所以即时问题与numpy。您使用了无效的 sympy 表达式！

===

您的basis函数显示相同问题的一个方面。我们需要再次查看 FULL 回溯，然后测试部分以确定确切的问题表达方式。

In [16]: basis(x, -1,0.5,0)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-16-b328f95b3c79> in <module>
----> 1 basis(x, -1,0.5,0)

<ipython-input-15-c6436540e3f3> in basis(x, a, b, h)
      1 def basis(x,a,b, h):
----> 2     if x <= a-2*h:
      3         return 0;
      4     elif (x<=a-h) or (x >=2*h):
      5         return (1/6)*(2*h+(x-a))**3

/usr/local/lib/python3.8/dist-packages/sympy/core/relational.py in __nonzero__(self)
    382 
    383     def __nonzero__(self):
--> 384         raise TypeError("cannot determine truth value of Relational")
    385 
    386     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

此表达式是sympy关系：

In [17]: x <= -1
Out[17]: x ≤ -1

但我们不能在Python if语句中使用这种关系。

In [18]: if x <= -1: pass
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-18-b56148a48367> in <module>
----> 1 if x <= -1: pass

/usr/local/lib/python3.8/dist-packages/sympy/core/relational.py in __nonzero__(self)
    382 
    383     def __nonzero__(self):
--> 384         raise TypeError("cannot determine truth value of Relational")
    385 
    386     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

Python if 是简单的真/假开关；它的论点必须评估为一个或另一个。错误正在告诉我们，sympy.Relational不起作用。 0 < x < 1 是基本 Python if 的变体（它测试 0<x 和 x<1 并执行 and）。

这是一个变体，我们经常在numpy（和pandas）是：

In [20]: 0 < np.array([0,1,2])
Out[20]: array([False,  True,  True])

In [21]: 0 < np.array([0,1,2])<1
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-21-bc1039cec1fc> in <module>
----> 1 0 < np.array([0,1,2])<1

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

numpy 表达式有多个 True/False 值，不能用于需要简单 True/False 的 Python 表达式。

编辑

正确扩展两侧测试：

In [23]: expr = np.piecewise((0, x-a <=  -2*h),
    ...:  ((1/6)*(2*h+(x-a))**3, (-2*h<=x-a)&(x-a<=-h)),
    ...:  (2*h**3/3-0.5*(x-a)**2*(2*h+(x-a)), (-h<= x-a)&(x-a<= 0)),
    ...:  (2*(h**3/3)-0.5*(x-a)**2*(2*h+(x-a)), (0<=x-a)&(x-a<=2*h)),
    ...:  ((1/6)*(2*h-(x-a))**3, (h<=x-a)&(x-a<=2*h)), (0, x-a<=2*h))

In [24]: expr
Out[24]: 
array([-0.5*(x + 1)**2*(x + 2.0) + 0.0833333333333333,
       -0.5*(x + 1)**2*(x + 2.0) + 0.0833333333333333], dtype=object)

In [26]: p = lambdify((x,), expr)

x是expr唯一的Sympy符号@。

查看生成的函数：

In [27]: print(p.__doc__)
Created with lambdify. Signature:

func(x)

Expression:

[-0.5*(x + 1)**2*(x + 2.0) + 0.0833333333333333  -0.5*(x + 1)**2*(x + 2.0)...

Source code:

def _lambdifygenerated(x):
    return ([-0.5*(x + 1)**2*(x + 2.0) + 0.0833333333333333, -0.5*(x + 1)**2*(x + 2.0) + 0.0833333333333333])