Would this be of any interest to anyone (see and run attached). I have
written a piecewise function class (not like the one already in sympy)
where the piecewise function is defined on an ordered numerical grid (see
attached). The class in not finished but I can calculate the convolved
powers of gate functions and will extend the class to calculate the
convolution of two piecewise functions if sympy can integrate the member
functions of the piecewise function. My interest comes from knowing the
Fourier transform of a gate function is a sinc function. Thus the Fourier
transform of a gate repeatedly convolved with itself is the power of a sinc
function. This is relevant to finding rf modulation envelopes that
minimize sidebands. Please run the code an see if it of general interest
and I will complete the class and documentation and submit it as an
addition to sympy.
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/c9634965-089a-4456-b27b-ff3f0211a042n%40googlegroups.com.
from sympy import *
import bisect,copy
#init_printing( use_latex='mathjax' )
def sympy_fct(f,x,t):
return f.subs(x,t)
def union(a,b):
"""
If a and b are ordered lists return an ordered list containing all
of the elements of a and b with no duplicates
"""
return sorted(list(set(a)|set(b)))
def interval(x,x0):
"""
If x is a sorted list of lenght n then
x0 < x[0] -> 0
x[i] < x0 <= x[i+1] -> i+1
x0 >= x[-1] -> n
"""
if x0 < x[0]:
return 0
if x0 >= x[-1]:
return len(x)
return bisect.bisect_left(x,x0)
def Integrate(f,xlo,xhi):
x = Symbol('x',real=True)
#print('Integrate: ',f,xlo,xhi)
if f == 0 or f == S(0):
return 0
return integrate(f,(x,xlo,xhi))
def Subs(f,x,xp):
"""
Takes care of the case where f is a python constant and not a sympy
expression.
"""
#print('Subs: ',f,x,xp)
if isinstance(f,Expr):
return f.subs(x,xp)
return f
class PieceWiseFunction:
"""
A piecewise function is defined by a ordered list of numbers self.x
and a list of sympy expressions (functions) self.f of the symbol x.
If the length of self.x is n then the length of self.f is n+1. Then
the fuction is given by:
x < self.x[0] -> self.f0[0]
x < self.x[i] and x >= self.x[i+1] -> self.f[i+1]
x >= self.x[-1] -> self.f[-1]
The init options are to input self.f and self.x separately or to
input one list [self.f[0],self.x[0],self.f[1],...,self.x[-1],self.f[-1]]
x = Symbol('x',real=True)
pwf1 = PieceWiseFunction([exp(x),-1,x**2,1,exp(-x)])
pwf2 = PieceWiseFunction([exp(x),x**2,exp(-x)],[-1,1])
pwf1 and pwf2 define the same piewcewise fuction:
exp(x): x < -1
x**2: -1 <= x < 1
exp(-x): x >= 1
"""
def __init__(self,pf,x=[]):
self.x = []
self.f = []
if len(pf) > 0 and len(x) == 0:
for i in range(len(pf)):
if i % 2 == 0:
if isinstance(pf[i],Expr):
self.f.append(pf[i])
else:
self.f.append(pf[i]*S(1))
else:
self.x.append(pf[i])
elif len(x) > 0:
self.x = x
self.f = pf
def fct(self,t):
"""
Return the value of the PieceWiseFunction for value t where t is
a number not a symbol.
"""
if self.x == []:
return(S(0))
if t < self.x[0]:
if isinstance(self.f[0],Expr):
return self.f[0].subs(x,t)
else:
return self.f[0]
elif t > self.x[-1]:
if isinstance(self.f[-1],Expr):
return self.f[-1].subs(x,t)
else:
return self.f[-1]
else:
i = bisect.bisect_left(self.x,t)-1
if isinstance(self.f[i+1],Expr):
return self.f[i+1].subs(x,t)
else:
return self.f[i+1]
def fct_expr(self,t):
"""
Return the expression of the PieceWiseFunction for the interval
that t is on where t is a number not a symbol.
"""
if self.x == []:
return(S(0))
if t < self.x[0]:
return self.f[0]
elif t > self.x[-1]:
return self.f[-1]
else:
i = bisect.bisect_left(self.x,t)-1
return self.f[i+1]
def RemoveDuplcates(self):
"""
Simplify a PieceWiseFunction that has the same expression defined
on adjacent intervals.
"""
i = 0
while i < len(self.x):
if self.f[i] == self.f[i+1]:
self.x.pop(i)
self.f.pop(i+1)
else:
i +=1
print(self.x)
print(self.f)
return
def expand(self):
"""
Expand the expression defined on every interval of a
PieceWiseFunction.
"""
f = []
for fi in self.f:
f.append(expand(fi))
self.f = f
return PieceWiseFunction(f,self.x)
def __mul__(self,a):
"""
Multiply two PieceWiseFunctions together. a can be another
PieceWiseFunction, a sympy expression, or a number. If a is
a expression or a number it multiplies all the functions defining
the PieceWiseFunction.
"""
if isinstance(a,(Expr,int,float)):
F = []
for f in self.f:
F.append(a*f)
b = PieceWiseFunction(F,self.x)
return b
if isinstance(a,PieceWiseFunction):
x = union(self.x,a.x)
f1 = self.fct_expr(x[0]-1)
f2 = a.fct_expr(x[0]-1)
f = [f1*f2]
xnew = [x[0]]
for i in range(0,len(x)-1):
xmid = (x[i]+x[i+1])/2
f1 = self.fct_expr(xmid)
f2 = a.fct_expr(xmid)
f1f2 = f1*f2
f.append(f1f2)
xnew.append(x[i])
xnew.append(x[i+1])
f1 = self.fct_expr(x[-1]+1)
f2 = a.fct_expr(x[-1]+1)
f.append(f1*f2)
xnew.append(x[-1])
xnew = sorted(list(set(xnew)))
pwf = PieceWiseFunction(f,xnew)
pwf.RemoveDuplcates()
return pwf
def __add__(self,a):
"""
Add two PieceWiseFunctions together. a can be another
PieceWiseFunction, a sympy expression, or a number. If a is
a expression or a number it is added to all the functions defining
the PieceWiseFunction.
"""
if isinstance(a,(Expr,int,float)):
F = []
for f in self.f:
F.append(a+f)
b = PieceWiseFunction(self.x,F)
return b
if isinstance(a,PieceWiseFunction):
x = union(self.x,a.x)
f1 = self.fct_expr(x[0]-1)
f2 = a.fct_expr(x[0]-1)
f = [f1+f2]
xnew = [x[0]]
for i in range(0,len(x)-1):
xmid = (x[i]+x[i+1])/2
f1 = self.fct_expr(xmid)
f2 = a.fct_expr(xmid)
f1f2 = f1+f2
f.append(f1f2)
xnew.append(x[i])
xnew.append(x[i+1])
f1 = self.fct_expr(x[-1]+1)
f2 = a.fct_expr(x[-1]+1)
f.append(f1+f2)
xnew.append(x[-1])
xnew = sorted(list(set(xnew)))
pwf = PieceWiseFunction(f,xnew)
pwf.RemoveDuplcates()
return pwf
def __neg__(self):
negf = []
for f in self.f:
negf.append(-f)
return PieceWiseFunction(negf,self.x)
def __rmul__(self,a):
if isinstance(a,(Expr,int,float)):
return self*a
return self*a
def shift(self,t0):
"""
shift shifts the PieceWiseFunction t0 to the right. t0 must be
a number and not an expression.
"""
x = Symbol('x',real=True)
x0 = []
for xi in self.x:
x0.append(xi+t0)
f0 = []
for fi in self.f:
if isinstance(fi,Expr):
f0.append(fi.subs(x,x-t0).expand())
else:
f0.append(fi)
return PieceWiseFunction(f0,x0)
def __str__(self):
"""
Returns a string in list format where the string is
[{self.f[0]},self.x[0],...,{self.f[i]},self.x[i],...,self.x[-1],{self.f[-1]}]
and the braces {} distiguish fuctions from grid points.
"""
s = '[{'+str(self.f[0])+'},'
for i in range(1,len(self.f)):
s += str(self.x[i-1])+','
s += '{'+str(self.f[i])+'},'
s = s[:-1]+']'
return s
def _repr_latex_(self):
return str(self)
def __sub__(self,a):
return self + (-a)
def integral(self,xlims=[]):
"""
If xlims=[], integral returns a PieceWiseFunction that is the
integral of self between [-oo,x]. If xlims=[xlo,xhi] and xlo and
xhi are numbers integral returns the integral of self between
[xlo,xhi]. xlo and xhi must be numbers so integral can determine
which intervals of self.x they are in. A number or an expression
can be returned depending on the form of self.f.
"""
if xlims == []:
"""
Computes the running integral of a PieceWiseFunction
on [-oo,x] returning a PieceWiseFunction.
"""
x0 = Symbol('x0',real=True);
intf = integrate(self.f[0],(x,-oo,x0))
intf = intf.subs(x0,x)
F = [intf]
intf0 = intf.subs(x,self.x[0])
for i in range(1,len(self.x)):
intf = integrate(self.f[i],(x,self.x[i-1],x0))
intf = intf.subs(x0,x)+intf0
intf0 = intf.subs(x,self.x[i])
F.append(intf)
intf = integrate(self.f[-1],(x,self.x[-1],x0))
intf = intf.subs(x0,x)+intf0
F.append(intf)
return PieceWiseFunction(F,self.x)
else:
"""
Computes the integral of a PieceWiseFunction on the interval
[xlo,xhi], both xlo and xhi must be numbers not symbols or
expressions.
"""
xlo = xlims[0]
xhi = xlims[1]
nx = len(self.x)
x = Symbol('x',real=True)
"""
Intervals are numbered 0 if x <= self.x[0] and len(self.x)
if x > self.x[-1] and i if self.x[i-1] <= x < self.x[i].
"""
ixlo = interval(self.x,xlo) # self.x[] interval xlo is in
ixhi = interval(self.x,xhi) # self.x[] interval xhi is in
if ixlo == ixhi: # xlo and xhi in the same interval in self.x
S = integrate(self.f[ixlo],(x,xlo,xhi))
return S
S = 0
if ixlo == 0: # xlo <= self.x[0]
S += integrate(self.f[0],(x,xlo,self.x[0]))
if ixhi == len(self.x): # xhi > self.x[-1]
S += integrate(self.f[-1],(x,self.x[-1],xhi))
if ixlo == 0 and ixhi == nx:
for i in range(nx-1):
S += integrate(self.f[i+1],(x,self.x[i],self.x[i+1]))
return S
if ixlo > 0:
S += integrate(self.f[ixlo],(x,xlo,self.x[ixlo]))
if ixhi < nx:
S += integrate(self.f[ixhi],(x,self.x[ixhi-1],xhi))
if ixhi-ixlo > 1:
for i in range(ixlo+1,ixhi):
S += integrate(self.f[i],(x,self.x[i-1],self.x[i]))
return S
def remapintervals(self,x0):
"""
For each interval [self.x[i],self.x[i+1]] determines interval
self.x[i]-x0 snf self.x[i+1]-x0 is in (ixlo,ixhi). If both points
are not in the same interval map all intermediate grid points
back to [self.x[i],self.x[i+1]] by adding x0 and store values in
a list.
"""
nx = len(self.x)
X = []
for i in range(nx-1):
ixlo = interval(self.x,self.x[i]-x0)
ixhi = interval(self.x,self.x[i+1]-x0)
if ixhi > ixlo:
for j in range(ixlo,ixhi):
xtmp = self.x[j]+x0
if xtmp not in self.x:
X.append(xtmp)
return X
def subs(self,x,y):
F = []
for f in self.f:
if isinstance(f,Expr):
F.append(f.subs(x,y))
else:
F.append(f)
return PieceWiseFunction(F,self.x)
def refinegrid(self,X0):
"""
X0 is an ordered list of new grid points. Insert the grid points
in X0 to a copy of self.x maintaining the order in the copy and
insert in a copy of self.f in the position
same position
"""
X = self.x.copy()
F = self.f.copy()
for x0 in X0:
if x0 not in self.x:
ix0 = interval(X,x0)
f0 = F[ix0]
X = X[:ix0]+[x0]+X[ix0:]
F = F[:ix0]+[f0]+F[ix0:]
return PieceWiseFunction(F,X)
def gateintegral(self,L):
"""
Integrate PieceWiseFunction between limits [x-L,x] and
return a PieceWiseFunction (convolution of self with gate of
length L). L must be a number and not a sympy expression.
"""
(x,xt) = symbols(r'x \tilde{x}',real=True)
# Add self.x[i]+L to PieceWiseFunction grid
Xnew = []
for xi in self.x:
xpL = xi+L
if xpL not in self.x:
Xnew.append(xpL)
pwf = self.refinegrid(Xnew)
# Gate integral for x < pwf.x[0] on refined grid
f = Integrate(pwf.f[0],xt-L,xt)
X = [pwf.x[0]]
F = [f]
# Gate integral for pwf.x[0] <= x < pwf.x[-1] on refined grid
for i in range(1,len(pwf.x)):
j = interval(pwf.x,pwf.x[i]-L)
if j == i: # x and x-L in same interval
f = Integrate(pwf.f[i],xt-L,xt)
X.append(pwf.x[i])
F.append(f)
if j+1 == i: # x and x-L in adjacent intervals
f = Integrate(pwf.f[j],xt-L,pwf.x[j]) + \
Integrate(pwf.f[i],pwf.x[j],xt)
X.append(pwf.x[i])
F.append(f)
if j+1 < i: # x and x-L in non-adjacent intervals
f = Integrate(pwf.f[j],xt-L,pwf.x[j]) + \
Integrate(pwf.f[i],pwf.x[i-1],xt) + \
pwf.integral([pwf.x[j],pwf.x[i-1]])
X.append(pwf.x[i])
F.append(f)
# Gate integral for pwf.x[-1] <= x on refined grid
f = Integrate(pwf.f[-1],xt-L,xt)
F.append(f)
pwf = PieceWiseFunction(F,X)
return pwf.subs(xt,x)
def asy(self,name):
"""
Exports an PieceWiseFunction as function writen in Asymptote code
for the purpose of using the Asymptote software to plot the
function.
https://asymptote.sourceforge.io/
"""
sp = ' '
i = 1
F = 'REALFCT[] F'+name+' = {'
X = 'real[] X'+name+' = '+str(self.x).replace('[','{').replace(']','}')+';\n'
s = '//Asymptote code for piecewise function '+name+'\n'
s+= '//generated by piecewise.py for data\n'
s+= '//'+str(self)+'\n'
s+= X
for f in self.f:
col = 0
s += 'real F'+name+str(i)+'(real x)\n{\n'
col += 2
s+= col*sp+'return '+str(f) +';\n'
col -= 2
s+= col*sp+'}\n'
F += 'rfct(F'+name+str(i)+'),'
i += 1
F = F[:-1]+'};\n'
return s+F+'PIECEWISE '+name+' = piecewise(X'+name+',F'+name+');\n'
"""
Test to calculate convolved powers of gate function h0. Look at doc
string for __str__ to understand output.
"""
h0 = PieceWiseFunction([S(0),S(1),S(0)],[-S(1)/S(2),S(1)/S(2)])
print('h0 = ',h0)
print('h0.integral([-oo,oo]) = ',h0.integral([-oo,oo]))
print('FourierTransform{h0} = sinc(f)')
h1 = h0.gateintegral(1)
h1 = h1.shift(-S(1)/S(2))
print('h1 = ',h1)
print('h1.integral([-oo,oo]) = ',h1.integral([-oo,oo]))
print('FourierTransform{h1} = sinc(f)**2')
h2 = h1.gateintegral(1)
h2 = h2.shift(-S(1)/S(2)).expand() #Center h2
print('h2 = ',h2)
print('h2.integral([-oo,oo]) = ',h2.integral([-oo,oo]))
print('FourierTransform{h2} = sinc(f)**3')
h3 = h2.gateintegral(1)
h3 = h3.shift(-S(1)/S(2)) #Center h3
print('h3 = ',h3)
print('h3.integral([-oo,oo]) = ',h3.integral([-oo,oo]))
print('FourierTransform{h3} = sinc(f)**4')
h4 = h3.gateintegral(1)
h4 = h4.shift(-S(1)/S(2)).expand() #Center h4
print('h4 = ',h4)
print('h4.integral([-oo,oo]) = ',h4.integral([-oo,oo]))
print('FourierTransform{h4} = sinc(f)**5')
