CMSIS-DSP/PythonWrapper/examples/testdsp5.py

import cmsisdsp as dsp 
import numpy as np
from numpy.testing import assert_allclose
from scipy.stats import entropy,tstd, tvar
from scipy.special import logsumexp
from scipy.linalg import cholesky,ldl,solve_triangular
from scipy import signal

def imToReal1D(a):
    ar=np.zeros(np.array(a.shape) * 2)
    ar[0::2]=a.real
    ar[1::2]=a.imag
    return(ar)

def realToIm1D(ar):
    return(ar[0::2] + 1j * ar[1::2])

print("Max and AbsMax")
a=np.array([1.,-3.,4.,0.,-10.,8.])

i=dsp.arm_absmax_no_idx_f32(a)
print(i)
assert i==10.0

i=dsp.arm_absmax_no_idx_f64(a)
print(i)
assert i==10.0

r,i=dsp.arm_absmax_f64(a)
assert i==4
assert r==10.0

r,i=dsp.arm_max_f64(a)
assert i==5
assert r==8.0

i=dsp.arm_max_no_idx_f32(a)
print(i)
assert i==8

i=dsp.arm_max_no_idx_f64(a)
print(i)
assert i==8

print("Min and AbsMin")
a=np.array([1.,-3.,4.,0.5,-10.,8.])

i=dsp.arm_absmin_no_idx_f32(a)
print(i)
assert i==0.5

i=dsp.arm_absmin_no_idx_f64(a)
print(i)
assert i==0.5

r,i=dsp.arm_absmin_f64(a)
assert i==3
assert r==0.5

r,i=dsp.arm_min_f64(a)
assert i==4
assert r==-10

i=dsp.arm_min_no_idx_f32(a)
print(i)
assert i==-10

i=dsp.arm_min_no_idx_f64(a)
print(i)
assert i==-10

print("Barycenter")

a=np.zeros((12,3))
w=np.array([[2.0] * 12])[0]
w[0]=3
w[11]=3
a[0] =[0., 0., -0.951057]
a[1] =[0., 0., 0.951057]
a[2] =[-0.850651, 0., -0.425325]
a[3] =[0.850651, 0., 0.425325]
a[4] =[0.688191, -0.5, -0.425325]
a[5] =[0.688191, 0.5, -0.425325]
a[6] =[-0.688191, -0.5, 0.425325]
a[7] =[-0.688191, 0.5, 0.425325]
a[8] =[-0.262866, -0.809017, -0.425325]
a[9] =[-0.262866, 0.809017, -0.425325]
a[10]=[0.262866, -0.809017, 0.425325]
a[11]=[0.262866, 0.809017, 0.425325]


scaled= np.dot(a.T , w)
print(scaled)
ref=scaled/np.sum(w)
#print(ref)

points = np.array(a).reshape(12*3)
weights = w.reshape(12)

#print(points)
#print(weights)

result=dsp.arm_barycenter_f32(points,weights,12,3)


assert_allclose(ref,result,rtol=1e-6,atol=1e-6)

print("Weighted sum")

nb=10
s = np.random.randn(nb)
w = np.random.randn(nb)

ref=np.dot(s,w)/np.sum(w)
print(ref)

res=dsp.arm_weighted_sum_f32(s,w)
print(res)

assert_allclose(ref,res,2e-5)

print("Entropy")
s = np.abs(np.random.randn(nb))
s = s / np.sum(s)

ref=entropy(s)
print(ref)
res=dsp.arm_entropy_f32(s)
print(res)
assert_allclose(ref,res,1e-6)

res=dsp.arm_entropy_f64(s)
print(res)
assert_allclose(ref,res,1e-10)

print("Kullback-Leibler")
sa = np.abs(np.random.randn(nb))
sa = sa / np.sum(sa)

sb = np.abs(np.random.randn(nb))
sb = sb / np.sum(sb)

ref=entropy(sa,sb)
print(ref)
res=dsp.arm_kullback_leibler_f32(sa,sb)
print(res)
assert_allclose(ref,res,1e-6)

res=dsp.arm_kullback_leibler_f64(sa,sb)
print(res)
assert_allclose(ref,res,1e-10)

print("Logsumexp")
s = np.abs(np.random.randn(nb))
s = s / np.sum(s)

ref=logsumexp(s)
print(ref)
res=dsp.arm_logsumexp_f32(s)
print(res)
assert_allclose(ref,res,1e-6)

print("Logsumexp dot prod")
sa = np.abs(np.random.randn(nb))
sa = sa / np.sum(sa)

sb = np.abs(np.random.randn(nb))
sb = sb / np.sum(sb)

d = 0.001
# It is a proba so must be in [0,1]
# But restricted to ]d,1] so that the log exists
sa = (1-d)*sa + d
sb = (1-d)*sb + d

ref=np.log(np.dot(sa,sb))
print(ref)

sa = np.log(sa)
sb = np.log(sb)

res=dsp.arm_logsumexp_dot_prod_f32(sa,sb)
print(res)
assert_allclose(ref,res,3e-6)

print("vexp")
sa = np.random.randn(nb)

ref = np.exp(sa)
print(ref)

res=dsp.arm_vexp_f32(sa)
print(res)
assert_allclose(ref,res,1e-6)

res=dsp.arm_vexp_f64(sa)
print(res)
assert_allclose(ref,res,1e-10)


print("vlog")
sa = np.abs(np.random.randn(nb)) + 0.001

ref = np.log(sa)
print(ref)

res=dsp.arm_vlog_f32(sa)
print(res)
assert_allclose(ref,res,2e-5,1e-5)

res=dsp.arm_vlog_f64(sa)
print(res)
assert_allclose(ref,res,2e-9,1e-9)

print("Cholesky")

a=np.array([[4,12,-16],[12,37,-43],[-16,-43,98]])
ref=cholesky(a,lower=True)
print(ref)

status,res=dsp.arm_mat_cholesky_f32(a)
print(res)
assert_allclose(ref,res,1e-6,1e-6)

status,res=dsp.arm_mat_cholesky_f64(a)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("LDLT")

def swaprow(m,k,j):
    tmp = np.copy(m[j,:])
    m[j,:] = np.copy(m[k,:])
    m[k,:] = tmp 
    return(m)

# F32 test
status,resl,resd,resperm=dsp.arm_mat_ldlt_f32(a)
n=3
p=np.identity(n)
for k in range(0,n):
    p = swaprow(p,k,resperm[k])

res=resl.dot(resd).dot(resl.T)

permutedSrc=p.dot(a).dot(p.T)
print(res)
print(permutedSrc)

assert_allclose(permutedSrc,res,1e-5,1e-5)

# F64 test
print("LDLT F64")
status,resl,resd,resperm=dsp.arm_mat_ldlt_f64(a)
n=3
p=np.identity(n)
for k in range(0,n):
    p = swaprow(p,k,resperm[k])

res=resl.dot(resd).dot(resl.T)

permutedSrc=p.dot(a).dot(p.T)
print(res)
print(permutedSrc)

assert_allclose(permutedSrc,res,1e-9,1e-9)

print("Solve lower triangular")
a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 1, 1, 1]])
b = np.array([[4,2,4,2],[8,4,8,4]]).T
x = solve_triangular(a, b,lower=True)
print(a)
print(b)
print(x)

b = np.array([[4,2,4,2],[8,4,8,4]]).T
status,res=dsp.arm_mat_solve_lower_triangular_f32(a,b)
print(res)
assert_allclose(x,res,1e-5,1e-5)


b = np.array([[4,2,4,2],[8,4,8,4]]).T
status,res=dsp.arm_mat_solve_lower_triangular_f64(a,b)
print(res)
assert_allclose(x,res,1e-9,1e-9)


print("Solve upper triangular")
a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 1, 1, 1]])
b = np.array([[4,2,4,2],[8,4,8,4]]).T
x = solve_triangular(a.T, b,lower=False)
print(a.T)
print(b)
print(x)

b = np.array([[4,2,4,2],[8,4,8,4]]).T
status,res=dsp.arm_mat_solve_upper_triangular_f32(a.T,b)
print(res)
assert_allclose(x,res,1e-5,1e-5)


b = np.array([[4,2,4,2],[8,4,8,4]]).T
status,res=dsp.arm_mat_solve_upper_triangular_f64(a.T,b)
print(res)
assert_allclose(x,res,1e-9,1e-9)


print("Mat mult f64")
a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 1, 1, 1]])
b = np.array([[4,2,4,2],[8,4,8,4]]).T 

ref =a.dot(b)
print(ref)

status,res = dsp.arm_mat_mult_f64(a,b)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("mat sub f64")
a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 1, 1, 1]])
b = a.T 

ref = a - b
print(ref)

status,res = dsp.arm_mat_sub_f64(a,b)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("abs f64")
s = np.random.randn(nb)
ref = np.abs(s)
res=dsp.arm_abs_f64(s)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("add f64")
sa = np.random.randn(nb)
sb = np.random.randn(nb)
ref = sa + sb
res=dsp.arm_add_f64(sa,sb)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("sub f64")
sa = np.random.randn(nb)
sb = np.random.randn(nb)
ref = sa - sb
res=dsp.arm_sub_f64(sa,sb)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("dot prod  f64")
sa = np.random.randn(nb)
sb = np.random.randn(nb)
ref = sa.dot(sb)
res=dsp.arm_dot_prod_f64(sa,sb)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("mult f64")
sa = np.random.randn(nb)
sb = np.random.randn(nb)
ref = sa * sb
res=dsp.arm_mult_f64(sa,sb)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("negate f64")
sa = np.random.randn(nb)
ref = -sa
res=dsp.arm_negate_f64(sa)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("offset f64")
sa = np.random.randn(nb)
ref = sa + 0.1
res=dsp.arm_offset_f64(sa,0.1)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("scale f64")
sa = np.random.randn(nb)
ref = sa * 0.1
res=dsp.arm_scale_f64(sa,0.1)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("mean f64")
sa = np.random.randn(nb)
ref = np.mean(sa)
res=dsp.arm_mean_f64(sa)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("power f64")
sa = np.random.randn(nb)
ref = np.sum(sa * sa)
res=dsp.arm_power_f64(sa)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("std f64")
sa = np.random.randn(nb)
ref = tstd(sa)
res=dsp.arm_std_f64(sa)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("variance f64")
sa = np.random.randn(nb)
ref = tvar(sa)
res=dsp.arm_var_f64(sa)
print(ref)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("fill f64")
nb=20
ref = np.ones(nb)*4.0 
res = dsp.arm_fill_f64(4.0,nb)
assert_allclose(ref,res,1e-10,1e-10)

print("copy f64")
nb=20
sa = np.random.randn(nb)
ref = sa
res = dsp.arm_copy_f64(sa)
assert_allclose(ref,res,1e-10,1e-10)

print("arm_div_int64_to_int32")
den=0x7FFF00000000 
num=0x10000
ref=den//num

res=dsp.arm_div_int64_to_int32(den,num)
print(ref)
print(res)

print("fir f64")

firf64 = dsp.arm_fir_instance_f64()
dsp.arm_fir_init_f64(firf64,3,[1.,2,3],[0,0,0,0,0,0,0])
filtered_x = signal.lfilter([3,2,1.], 1.0, [1,2,3,4,5,1,2,3,4,5])
print(filtered_x)
ra=dsp.arm_fir_f64(firf64,[1,2,3,4,5])
rb=dsp.arm_fir_f64(firf64,[1,2,3,4,5])
assert ((filtered_x == np.hstack([ra,rb])).all)

print("arm_cmplx_mag")
sa = np.random.randn(nb)
ca = realToIm1D(sa)
ref = np.abs(ca)
print(ref)
res=dsp.arm_cmplx_mag_f32(sa)
print(res)
assert_allclose(ref,res,1e-6,1e-6)
res=dsp.arm_cmplx_mag_f64(sa)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("arm_cmplx_mag_squared")
sa = np.random.randn(nb)
ca = realToIm1D(sa)
ref = np.abs(ca) * np.abs(ca)
print(ref)
res=dsp.arm_cmplx_mag_squared_f32(sa)
print(res)
assert_allclose(ref,res,1e-6,1e-6)
res=dsp.arm_cmplx_mag_squared_f64(sa)
print(res)
assert_allclose(ref,res,1e-10,1e-10)

print("cmplx mult")
sa = np.random.randn(nb)
ca = realToIm1D(sa)
sb = np.random.randn(nb)
cb = realToIm1D(sb)
ref = imToReal1D(ca * cb)
print(ref)
res = dsp.arm_cmplx_mult_cmplx_f32(sa,sb)
print(res)
assert_allclose(ref,res,1e-6,1e-6)

res = dsp.arm_cmplx_mult_cmplx_f64(sa,sb)
print(res)
assert_allclose(ref,res,1e-10,1e-10)