2bit_tests_git/fits_stats.py at master · tcromartie/2bit_tests_git · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
import numpy as np
import ctypes
import matplotlib.pyplot as plt
import scipy.stats
from astropy.io import fits

hdu = fits.open('guppi_55179_GBNCC06605_0045_short.fits',memmap=True)
data = hdu[1].data
z = np.squeeze(data[0]['DATA']).T
#print 'Array shape = %f' %(z.shape)

m = ctypes.c_double()
v = ctypes.c_double()
s = ctypes.c_double()
k = ctypes.c_double()

stats = ctypes.cdll.LoadLibrary("./stats.so")

def chan_stats(x,m,v,s,k):
    xptr = x.ctypes.data_as(ctypes.POINTER(ctypes.c_float))
    stats_out = stats.stats(xptr, len(z), ctypes.byref(m),ctypes.byref(v), ctypes.byref(s), ctypes.byref(k))
    #print 'mean = %f, variance = %f, skew = %f, kurtosis = %f' %(m.value, v.value, s.value, k.value)
    #print 'expected values: %f, %f, %f, %f' %(mu, sig, scipy.stats.skew(x), scipy.stats.kurtosis(x))
    #norm_kurt = scipy.stats.kurtosistest(x)
    #print norm_kurt
    #count, bins, ignored = plt.hist(x, 500, normed=True)
    #plt.plot(bins, 1/(sig * np.sqrt(2 * np.pi)) * np.exp(-(bins-mu)**2 / (2*sig**2)))
    #plt.show()
    return np.array([m.value, v.value, s.value, k.value])

istats = np.zeros((len(z),4))

for i in range(0,len(z)):
    x = z[i,:]
    istats[i] = chan_stats(x,m,v,s,k)

print istats
print istats[2000]