In [1]:
import tensorflow as tf
import numpy
import scipy.io
import glob
from DoppCorr import doppcorr
from GPUHough import frequencyHough
from CandSel import select_candidates
#first of all i define some job parameters:
#FFT length, security belt
tFft = 8192
#time interval of data
tObs = 9 #months
tObs = tObs*30*24*60*60 #seconds
#number of frequency belts in the hough where candidates are selected
numCands = 100
#security belt to avoid negative numbers in frequency-hough transform
securbelt = 4000
#frequency steps
stepFreq = 1/tFft
enhancement = 10
refinedStepFreq = stepFreq/enhancement
#epoch, defined as average time fo the observation run
epoch = (57722+57990)/2
#spindown range
spindownMin = -1e-9
spindownMax = 1e-10
stepSpindown = stepFreq/tObs
numStepsSpindown = numpy.round((spindownMax-spindownMin)/stepSpindown).astype(numpy.int64)
#and spindowns values
spindowns = numpy.arange(0, numStepsSpindown)
spindowns = numpy.multiply(spindowns,stepSpindown)
spindowns = numpy.add(spindowns, spindownMin)
#Loading data
#file paths
#pathSquare = ("PATH/TO/square%d.mat" % tFft)
#pathPatch = ("PATH/TO/square%dEcl.mat" % tFft)
#pathFolder = 'PATH/TO/DATAFOLDER'
pathSquare = "/home/protoss/TESI/wn100bkp/quadHWI52.mat"
pathPatch = "/home/protoss/TESI/wn100bkp/quadHWI52Ecl.mat"
pathFolder = '/home/protoss/TESI/wn100bkp/data/*.mat'
#data for doppler correction
#square = scipy.io.loadmat(pathSquare)['square'].astype(numpy.float64)
#numPoints = square.shape[0]
#patch = scipy.io.loadmat(pathPatch)['squareEcl'].astype(numpy.float64)
square = scipy.io.loadmat(pathSquare)['quad'].astype(numpy.float64)
numPoints = square.shape[0]
patch = scipy.io.loadmat(pathPatch)['quadratoEclNew'].astype(numpy.float64)
#list of input files (with peakmaps)
fileslist = glob.glob(pathFolder)
fileslist = numpy.array(fileslist)
fileslist = numpy.sort(fileslist)
numFiles = fileslist.size
#initialization of candidates array
allCands = numpy.zeros((numFiles,numPoints,numCands*2,9))
#WARNING at the moment the hough transform is serialized between different input files
#AND between sky positions
#TODO fully parallelize
#for loop into every input file in the frequency band chosen (tFFT 8192/4096 s)
for ithFile in numpy.arange(numFiles):
#loading the ith .mat input file and the useful data
structure = scipy.io.loadmat(fileslist[ithFile])['job_pack_0']
times = structure['peaks'][0,0][0]#.astype(numpy.float32)
frequencies = structure['peaks'][0,0][1]#.astype(numpy.float32)
weights = (structure['peaks'][0,0][4]+1)#.astype(numpy.float32)
#per doppler corr
velocities = structure['basic_info'][0,0]['velpos'][0,0][0:3,:].astype(numpy.float64)
numTimes = structure['basic_info'][0,0]['ntim'][0,0][0,0]
indices = structure['basic_info'][0,0]['index'][0,0][0]
#for loop into every sky positions in the area chosen ([+-10 gal lon, +-10 gal lat])
#the full sky grid is calculated with SNAG function pss_optmap,
#converted in galactic coordinates, with SNAG function astro_coord_fix
#filtered with a simple python script
#reconverted in ecliptic coordinates
#converted in a sky rectangle wth SNAG function astro2rect
for point in numpy.arange(0,numPoints):
#WARNING i didn't define a complete tensorflow input pipeline (this is a quite difficult problem and in fast developement), so i used only constants and i'm forced to reset the graph and the tf.Session at every iteration of the loop
#TODO build a complete input pipeline, using only variables with the (unavoidable) iterations
session = tf.Session()
#doppler correction, simple numpy porting of the SNAG function hfdf_patch
#TODO develop a fully vectorial GPU version of doppler correction
pointRectangle = square[point]
indicesOpt = indices-1
freqInCorr,freqCorr, freqForHough = doppcorr(pointRectangle, indicesOpt, velocities)
numStepsFreq = numpy.ceil(securbelt+(numpy.amax(freqCorr)-numpy.amin(freqCorr) + stepFreq + 2*refinedStepFreq)/refinedStepFreq)
#now let's use TensorFlow
#defining tf constants to the hough transform
#WARNING i'm using 64bit data for times because 9 months with a step of 8192s exceed 32bit precision, this forces to use 64bit spindowns to run tf.matmul (tensorflow needs same data type tensors). This slows a bit the code with a GPU capable to run 64bit calculations (eg Tesla series), but can be a problem on other series (eg GeForce)
#TODO remove any 64bit data
timesTF = tf.constant(times,dtype=tf.float32)
weightsTF = tf.constant(weights,dtype=tf.float32)
spindownsTF = tf.constant(spindowns, dtype=tf.float32)
timesHM = timesTF-epoch
timesHM = ((timesHM)*3600*24/refinedStepFreq)
timesHM = tf.cast(timesHM, tf.float32)
weightsHM = tf.reshape(weightsTF,(1,tf.size(weightsTF)))
weightsHM = weightsHM[0]
freqTF = tf.constant(freqForHough, dtype = tf.float32)
#calculating the Hough transofrm, with a fully vectorial tensorflow function with GPU parallelization
nRows = tf.constant(numStepsSpindown, dtype=tf.int32)
nColumns = tf.cast(numStepsFreq, dtype=tf.int32)
houghmap = frequencyHough(freqTF,timesHM,weightsHM,spindownsTF)
hough = session.run(houghmap)
#candidates selection, simple numpy porting of the SNAG function hfdf_peak
#TODO develop a fully vectorial GPU version of candidates selection
allCands[ithFile,point] = select_candidates(numCands, freqInCorr, hough, patch[point])
tf.reset_default_graph()
session.close()
#rearranging the candidates selected to write them in a .mat file compatible with the
#SNAG functions for coincidence selections and follow-up
numAllCands = numpy.int64(allCands.size/9)
allCands = allCands.reshape(numAllCands,9)
nonzeroCands = numpy.nonzero(allCands[:,0])
finallCands = allCands[nonzeroCands]
finallCands = numpy.transpose(finallCands)
#pathOut = 'PATH/TO/candidates.mat'
pathOut = 'candidates.mat'
scipy.io.savemat(pathOut,{'cand' : finallCands})
In [2]:
from matplotlib import pyplot
%matplotlib qt
pyplot.imshow(hough, aspect=20)
Out[2]:
In [ ]: