SymPy - Kronecker Delta Function Evaluation

I am using SymPy for the numerical analysis of large sets of equations. Part of my equation contains a Kronecker Delta function acting as an impulse such that when q = 0 -> dirac_delta = 1, otherwise dirac_delta = 0. I need to perform this for values of q = - 10 -> +10 in integer steps of 1.

A simplified example of my code is:

import sympy as sp
import numpy as np

modules = ["numpy", "sympy"]

# Create symbols
q = sp.symbols('q', integer=True)

# Create P_q symbol as a function of q
P_q = sp.symbols('P_q', cls=sp.Function)
P_q = P_q(q)

# Define Equation
# Simplified example KroneckerDelta - when q = 0, P_q = 1, otherwise P_q = 0
P_q_eq = sp.Eq(P_q, sp.KroneckerDelta(0,q))
P_q = sp.KroneckerDelta(0,q)
display(P_q_eq)

# Create a lambda function for fast numerical calculation 
lam_P_q = sp.lambdify(q, P_q, modules)

# Define the values of q
num_points = 21
data = np.linspace(-10, 10, num_points, dtype=int)
#print(data)

ans = lam_P_q(data)
print(ans)

On run I receive an error:

ValueError Traceback (most recent call last) in 36 #print(data) 37 ---> 38 ans = lam_P_q(data) 39 print(ans)

in _lambdifygenerated(q) 1 def _lambdifygenerated(q): ----> 2 return ((1 if 0 == q else 0))

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

My understanding is that I need to add the .any() or .all() because the lambdify is comparing the array q, to a single value 0. So when I modify the input data with .any() or .all() it then returns a single value.

However, I require a response of 0 or 1 for each value of q - such that it is an impulse response depending on the value of q.

print(q)

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

I have tried to provide the "0" comparison value as an array of zeroes with equal length to q, but this did not work.

I know that this can be generated with the scipy signal function "signal.unit_impulse(21, 'mid')", but I am unsure how to implement this in the SymPy format for lambdify to output q as outlined just above. I have tried to create a custom module to replace the sp.KroneckerDelta function to perform this but couldn't get a working solution (likely due to a poor implementation):

def impulse(num, p):
    val = signal.unit_impulse(num, 'mid')
    j = p + (num-1)/2
    kron_p = int(val[j])
        
    return kron_p

kronecker = {'KroneckerDelta': impulse(21,q)}
modules = [kronecker, "numpy", "sympy"]

Do I need to substitute the values of q into the lambdify function instead - if so, how do I specify a range of values to substitute?

I feel like I am doing something fundamentally wrong in my approach and would appreciate help getting what I thought would be a relatively simple thing to do in SymPy to work. (I am quite new to SymPy and definitely still trying to get my head around it). Thanks.

Solution

This is a bug in lambdify. Please open a sympy issue:

https://github.com/sympy/sympy/issues

You can work around it by rewriting the KroneckerDelta to a Piecewise:

In [12]: P_q
Out[12]: 
δ   
 0,q

In [13]: P_q.rewrite(Piecewise)
Out[13]: 
⎧0  for q ≠ 0
⎨            
⎩1  otherwise

In [14]: f = lambdify(q, P_q.rewrite(Piecewise), modules)

In [15]: f(data)
Out[15]: 
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0.])

This is the generated code:

In [16]: import inspect

In [18]: print(inspect.getsource(f))
def _lambdifygenerated(q):
    return select([not_equal(q, 0),True], [0,1], default=nan)