BUG: Spectral kurtosis GPU not matching CPU gold version

[telegraphic]

Kurtosis kernel giving different answer to CPU version! Would like to understand why

import numpy as np
from astropy import units as u
from hyperseti.hyperseti import spectral_kurtosis

def sk_cpu(x, N=1):
    x_sum = x.sum(axis=0)
    x2_sum = (x**2).sum(axis=0)
    n = x.shape[0]
    return (N*n+1) / (n-1) * (n*(x2_sum / (x_sum*x_sum)) - 1)

# zero drift test, no normalization
noise = np.random.random(size=32*1*512).reshape((32,1,512)) * 2
test_data = np.zeros(shape=(32, 1, 512), dtype='float32')
test_data[:, 0, 255] = 100
test_data[:, 0, 100] = 100*np.sin(4*np.arange(0,32))
test_data += noise**2

metadata = {'fch1': 1000*u.MHz, 
        'dt': 1.0*u.s, 
        'df': 1.0*u.Hz}

sk = spectral_kurtosis(test_data, metadata)
skc = sk_cpu(test_data.squeeze(), N=2)

plt.subplot(2,1,1)
plt.plot(sk)
plt.plot(skc)
plt.subplot(2,1,2)
plt.plot(sk - skc)

[texadactyl]

issue_17.zip

Stdout: diff max=0.031, mean=0.0000617, std=0.0013797

When diff is replaced with err = np.abs((sk.get() - skc) / skc), stdout: err max=0.004, mean=0.0000079, std=0.0001693

When I ported spectral_kurtosis from hyperseti to CPU numpy code and re-ran, I get: diff max=0.000, mean=0.0000000, std=0.0000000 and the diff plot (3rd) is flat.

So, its not the spectral_kurtosis algorithm that is causing the differences. Something seems to be different in accuracy between cupy and numpy.

I do not believe that dedoppler_kurtosis_kernel is being invoked in your sample program unless I am missing something.

[david-macmahon]

This is probably not related to the problem at hand, but I'm a little confused. This definition of sk_cpu looks like its input is power and its computing kurtosis from power and power squared (I'm unclear on the extra N, n, and 1 parts, but that's another question):

 def sk_cpu(x, N=1):
     x_sum = x.sum(axis=0)
     x2_sum = (x**2).sum(axis=0)
     n = x.shape[0]
     return (N*n+1) / (n-1) * (n*(x2_sum / (x_sum*x_sum)) - 1)

 # zero drift test, no normalization
 noise = np.random.random(size=32*1*512).reshape((32,1,512)) * 2
 test_data = np.zeros(shape=(32, 1, 512), dtype='float32')
 test_data[:, 0, 255] = 100
 test_data[:, 0, 100] = 100*np.sin(4*np.arange(0,32))
 test_data += noise**2

Maybe I misunderstand how np.random works, but assuming it returns values in the range [0,1) it seems like test_data[:, 0, 100] will be bounded by -100 on the low side and +102 on the high side, which seems more like voltage data than power data. Is that valid input data for the sk_cpu function?

[telegraphic]

Closing as stale -- currently happy with SK performance