import pyfits
import pandas
import statsmodels
import mdp
import numpy as np
# pylab's already imported...
import agpy # need PCA tools
import aplpy
cubes.py requires pywcs for some subimage_integ,aper_wordl2pix,getspec, and coords_in_image
matplotlib.rcParams['figure.figsize'] = (12,7)
matplotlib.rcParams['savefig.dpi'] = 120
D = pyfits.getdata('/Users/adam/work/fastphot/TEST.fits')
RF = pyfits.getdata('/Users/adam/work/fastphot/Robberto_Found.fits')
D = D.astype('float')
DF = pandas.DataFrame(D)
cov = DF.cov()
imshow(cov)
<matplotlib.image.AxesImage at 0x106285ad0>
imshow(DF)
<matplotlib.image.AxesImage at 0x1062c0bd0>
# above shows that the first couple observations were wrong - maybe focusing? Maybe off-source? Probably off source, and just by chance overlapped enough to get included
D = D[:,30:] # exclude first 30 stars
# first observing block - correlations are more clear here...
firstchunk = D[18:516,:] # first dim is time
# filter out sections with too few data points
OKfirst = np.isnan(firstchunk).sum(axis=0) < 20
firstchunk_filtered = firstchunk[:,OKfirst]
FCF = pandas.DataFrame(firstchunk_filtered)
FCFcov = FCF.cov()
# filter the RA/Dec data (etc) als0
RFfirst = RF[30:][OKfirst]
imshow(log10(FCFcov))
<matplotlib.image.AxesImage at 0x1085f8050>
# PCA subtraction - remove 0'th component
efuncarr,covmat,evals,evects = agpy.PCA_tools.efuncs(np.nan_to_num(firstchunk_filtered),return_others=True)
efuncarr_blanked = efuncarr.copy()
efuncarr_blanked[:,0] = 0
firstchunk_subd = np.inner(efuncarr_blanked,evects)
figure(); imshow(log10(firstchunk_subd.real))
# there is a clear difference between first 300ish and the rest!
figure(); imshow(evects.real);
figure(); semilogy(evals.real)
[<matplotlib.lines.Line2D at 0x10ce1df50>]
This last plot is probably the most interesting. The first ~200 components are decently correlated, then then next 150 or so are effectively uncorrelated. Explored at the bottom: there is nothing of interest (unfortunately) in those last couple components. Since there's no reason to expect intrinsic phenomena to be correlated, any "real signal" would have to be in those components. Not sure what the dip from $10^{-11}$ to $10^{-27}$ around component 325 is, though.
plot(efuncarr[:,:2])
WARNING: ComplexWarning: Casting complex values to real discards the imaginary part
WARNING: ComplexWarning: Casting complex values to real discards the imaginary part [numpy.core.numeric]
[<matplotlib.lines.Line2D at 0x108d54450>, <matplotlib.lines.Line2D at 0x108d54e50>]
# find out which are most correlated with which component (presumably first 2
# account for most of what we see above, in the PCA-sub'd waterfall plot)
plot(evects[:,0])
[<matplotlib.lines.Line2D at 0x108d40a90>]
# select based on correlation with evect 1
set1 = evects[:,0] > -0.05
set1.shape
(359,)
# redo above with set 1
efuncarr2,covmat2,evals2,evects2 = agpy.PCA_tools.efuncs(np.nan_to_num(firstchunk_filtered[:,set1]),return_others=True)
efuncarr_blanked2 = efuncarr2.copy()
efuncarr_blanked2[:,0] = 0
firstchunk_set1_subd = np.inner(efuncarr_blanked2,evects2)
print firstchunk_set1_subd.shape,firstchunk_filtered.shape
(498, 212) (498, 359)
figure(figsize=[12,7],dpi=300)
imshow(log10(firstchunk_set1_subd.real+1000),aspect=1/3.,vmin=2.5)
<matplotlib.image.AxesImage at 0x10cc7f510>
plot(firstchunk_set1_subd[:,100])
[<matplotlib.lines.Line2D at 0x110495410>]
subplot(311); plot(efuncarr2[:,0])
subplot(312); plot(efuncarr2[:,1])
subplot(313);
plot(firstchunk_set1_subd[:,47]) # this one probably belongs in the other set
[<matplotlib.lines.Line2D at 0x1104adc90>]
plot(firstchunk_set1_subd[:,47:51])
[<matplotlib.lines.Line2D at 0x1126c7250>, <matplotlib.lines.Line2D at 0x1126c73d0>, <matplotlib.lines.Line2D at 0x1126c7550>, <matplotlib.lines.Line2D at 0x1126c76d0>]
# redo above with anti-set 1
efuncarr3,covmat3,evals3,evects3 = agpy.PCA_tools.efuncs(np.nan_to_num(firstchunk_filtered[:,True-set1]),return_others=True)
efuncarr_blanked3 = efuncarr3.copy()
efuncarr_blanked3[:,0] = 0
firstchunk_set2_subd = np.inner(efuncarr_blanked3,evects3)
print firstchunk_set2_subd.shape,firstchunk_filtered.shape
(498, 147) (498, 359)
imshow(log10(firstchunk_set2_subd.real+1000),aspect=1/5.) # add 1000 so it shows up...
figure()
plot(firstchunk_set2_subd[:,5])
[<matplotlib.lines.Line2D at 0x1129aded0>]
# the subd image shows another subset/breakdown...
plot(efuncarr2[:,:5]*-1)
figure()
plot(efuncarr3[:,:5])
[<matplotlib.lines.Line2D at 0x115a48350>, <matplotlib.lines.Line2D at 0x115a484d0>, <matplotlib.lines.Line2D at 0x115a48650>, <matplotlib.lines.Line2D at 0x115a487d0>, <matplotlib.lines.Line2D at 0x115a48950>]
import warnings
warnings.filterwarnings('ignore')
import pywcs
FF = aplpy.FITSFigure('/Users/adam/work/orion/orion_finder_mosaic.fits')
FF.show_grayscale(pmin=5,pmax=95)
FF.show_markers(RFfirst[set1]['RAJ2000'],RFfirst[set1]['DEJ2000'])
FF.show_markers(RFfirst[True-set1]['RAJ2000'],RFfirst[True-set1]['DEJ2000'],edgecolor='b')
INFO: Auto-setting vmin to 1.537e+04 INFO: Auto-setting vmax to 1.716e+04
# What if we just get rid of ALL of the correlated components?
# PCA subtraction - remove first 200 components
efuncarr4,covmat4,evals4,evects4 = agpy.PCA_tools.efuncs(np.nan_to_num(firstchunk_filtered),return_others=True)
efuncarr_blanked4 = efuncarr4.copy()
efuncarr_blanked4[:,:200] = 0
firstchunk_subd_huge = np.inner(efuncarr_blanked4,evects4)
LS = np.log10(firstchunk_subd_huge.real+1000) #logscale, make all above zero
imshow(LS,vmin=np.percentile(LS,1),vmax=np.percentile(LS,99))
<matplotlib.image.AxesImage at 0x1138df4d0>
for ii in xrange(6):
subplot(6,1,ii+1)
plot(firstchunk_subd_huge[:,ii])
for ii in xrange(6):
subplot(6,1,ii+1)
plot(firstchunk_subd_huge[:,105+ii])
FTFC = np.fft.fft(firstchunk_subd_huge,axis=0)
for ii in xrange(6):
subplot(6,1,ii+1)
plot(np.fft.fftshift(FTFC[:,ii])[:FTFC.shape[0]/2])
for ii in xrange(6):
subplot(6,1,ii+1)
plot(np.fft.fftshift(FTFC[:,105+ii])[:FTFC.shape[0]/2])
for ii in xrange(6):
subplot(6,1,ii+1)
plot(firstchunk_subd_huge[:,105+ii])
plot(np.fft.fft(firstchunk_subd_huge[:,105+ii]))
LSFT = np.log10(FTFC.real+1000) #logscale, make all above zero
imshow(LSFT,vmin=np.percentile(LSFT,1),vmax=np.percentile(LSFT,99))
<matplotlib.image.AxesImage at 0x11233c650>
errorbar(arange(LSFT.shape[1]),LSFT.max(axis=0),yerr=LSFT.std(axis=0),linestyle='none')
<Container object of 3 artists>
plot(LSFT.max(axis=0)/LSFT.std(axis=0))
possibleperiod = LSFT.max(axis=0)/LSFT.std(axis=0) > 5e5
for ii in xrange(6):
subplot(6,1,ii+1)
plot(np.fft.fftshift(FTFC[:,possibleperiod][:,ii])[:FTFC.shape[0]/2])
# second observing block - correlations are more clear here...
secondchunk = D[516:1020,:] # second dim is time
# filter out sections with too few data points
OKsecond = np.isnan(secondchunk).sum(axis=0) < 10
secondchunk_filtered = secondchunk[:,OKsecond]
SCF = pandas.DataFrame(secondchunk_filtered)
SCFcov = SCF.cov()
# filter the RA/Dec data (etc) als0
RFsecond = RF[30:][OKsecond]
imshow(SCFcov)
<matplotlib.image.AxesImage at 0x1122747d0>
# redo above with anti-set 1
efuncarr5,covmat5,evals5,evects5 = agpy.PCA_tools.efuncs(np.nan_to_num(secondchunk_filtered),return_others=True)
efuncarr_blanked5 = efuncarr5.copy()
efuncarr_blanked5[:,0] = 0
secondchunk_subd = np.inner(efuncarr_blanked5,evects5)
print secondchunk_subd.shape,secondchunk_filtered.shape
(504, 215) (504, 215)
semilogy(evals5)
[<matplotlib.lines.Line2D at 0x116858d90>]
imshow(log10(secondchunk_subd.real+1000),aspect=1/3.)
<matplotlib.image.AxesImage at 0x112628290>
efuncarr_blanked5[:,:7] = 0
secondchunk_subd = np.inner(efuncarr_blanked5,evects5)
print secondchunk_subd.shape,secondchunk_filtered.shape
LS = (log10(secondchunk_subd.real+1000))
imshow(LS,vmin=np.percentile(LS,1),vmax=np.percentile(LS,99),aspect=1/3.)
(504, 215) (504, 215)
<matplotlib.image.AxesImage at 0x1123104d0>
print firstchunk_filtered.shape,efuncarr.shape,FCF.shape,FCF.mean(axis=0).shape
(498, 359) (498, 359) (498, 359) (359,)
atmonorm = efuncarr[:,0]/efuncarr[:,0].mean()
firstchunk_normed = firstchunk_filtered / array(FCF.mean(axis=0)) / atmonorm[:,newaxis]
nplots=10
for ii in xrange(nplots):
subplot(nplots,1,ii+1)
plot(firstchunk_normed[:,ii])
imshow(firstchunk_normed.real,vmin=0.5,vmax=2,aspect=0.5)
<matplotlib.image.AxesImage at 0x11dc49610>
spiky = where(firstchunk_normed.max(axis=0) > 2.5)[0]
dippy = where(firstchunk_normed.min(axis=0) < 0.1)[0]
print dippy
[208 298 300 301 304 310 311 320 323 331 332 338 340 341]
#help(scatter)
# lots of overplots
import glob
fitsfiles = glob.glob('/Users/adam/observations/tcam/UT121105/proc-t0[0-9]50.corr.fits')
FFlist = []
for fn in fitsfiles:
FF = aplpy.FITSFigure(fn)
FF.show_grayscale(pmin=1,pmax=99)
FF.show_markers(RFfirst[set1]['RAJ2000'],RFfirst[set1]['DEJ2000'])
FF.show_markers(RFfirst[True-set1]['RAJ2000'],RFfirst[True-set1]['DEJ2000'],edgecolor='b')
FF.show_markers(RFfirst[spiky]['RAJ2000'],RFfirst[spiky]['DEJ2000'],edgecolor=(0,1,0),marker='x')
FF.show_markers(RFfirst[dippy]['RAJ2000'],RFfirst[dippy]['DEJ2000'],edgecolor=(1,0,1),marker='x')
FFlist.append(FF)
INFO: Auto-setting vmin to 1.653e+04 INFO: Auto-setting vmax to 1.843e+04 INFO: Auto-setting vmin to 1.895e+04 INFO: Auto-setting vmax to 2.032e+04 INFO: Auto-setting vmin to 1.485e+04 INFO: Auto-setting vmax to 1.581e+04 INFO: Auto-setting vmin to 1.522e+04 INFO: Auto-setting vmax to 1.631e+04 INFO: Auto-setting vmin to 1.483e+04 INFO: Auto-setting vmax to 1.597e+04 INFO: Auto-setting vmin to 1.556e+04 INFO: Auto-setting vmax to 1.676e+04 INFO: Auto-setting vmin to 1.508e+04 INFO: Auto-setting vmax to 1.620e+04 INFO: Auto-setting vmin to 1.429e+04 INFO: Auto-setting vmax to 1.538e+04 INFO: Auto-setting vmin to 1.468e+04 INFO: Auto-setting vmax to 1.620e+04