In [ ]:

    

import tensorflow as tf
import numpy
import scipy.io
from tensorflow.python.client import timeline
import time
import glob


# PREPARO FUNZIONI

#TODO STUDIARE FATTIBILITÀ DI USARE TUTTO A 32BYTES

# calcola la correzione doppler per ogni punto del cielo
def doppcorr(i):
	quadratoNP = quadrato[i]
	indicesOpt = indices-1
	inizi = indicesOpt[:-1]
	fini = indicesOpt[1:]

	velocitas = numpy.zeros((3,frequenze.size))
	for i in numpy.arange(0,nTempi-1):
		velocitas[:,inizi[i]:fini[i]+1] = veloc[:,i:i+1]

	velPerPosIndex = numpy.dot(quadratoNP,velocitas)
	divisoreIndex = 1+velPerPosIndex
	freqCorr = frequenze / divisoreIndex
	##print(freqCorr)

	#mi ricavo l'header per le frequenze
	freqMin = numpy.amin(freqCorr)
	#freqMax = tf.reduce_max(freqCorr)
	freqIniz = freqMin- stepFreq/2 - stepFreqRaffinato

	freqFinal = freqCorr-freqIniz
	freqFinal = (freqFinal/stepFreqRaffinato)-round(enhancement/2+0.001)

	#freqs = tf.concat([[freqIniz], freqCorr], 0)

	return freqIniz, freqCorr, freqFinal#, nstepFrequenze

def noncorr():
	freqMin = numpy.amin(frequenze)
	freqIniz = freqMin - stepFreq/2 -stepFreqRaffinato
	freqNonCor = (frequenze -freqIniz)/stepFreqRaffinato-round(enhancement/2+0.001)
	freqNonCor = tf.constant(freqNonCor, dtype = tf.float64)
	return freqNonCor


# calcola la hough per ogni punto del cielo (per ogni spindown)
def inDaHough(i, freqHM):
	def houghizza(stepIesimo):
		sdTimed = tf.multiply(spindownsTF[stepIesimo], tempiHM, name = "Tdotpert")
		#sdTimed = tf.cast(sdTimed, dtype=tf.float32)

		appoggio = tf.round(freqHM-sdTimed+securbelt/2, name = "appoggioperindici")
		appoggio = tf.cast(appoggio, dtype=tf.int32)

		#pesi32 = tf.cast(pesiHM, dtype = tf.float32)

		valorisx = tf.unsorted_segment_sum(pesiHM, appoggio, nColumns)
		valorisx = tf.cast(valorisx, dtype=tf.float32)
		return valorisx

	houghDiff = tf.map_fn(houghizza, tf.range(0, nRows), dtype=tf.float32, parallel_iterations=8)

	def sliceInt():
		#faccio integrazione finale (vecchia versione senza conv)
		semiLarghezza = tf.round(enhancement/2+0.001)
		semiLarghezza = tf.cast(semiLarghezza, tf.int64)
		houghInt = houghDiff[:,enhancement:nColumns]-houghDiff[:,0:nColumns - enhancement]
		houghInt = tf.concat([houghDiff[:,0:enhancement],houghInt],1)
		return houghInt

	hough = sliceInt()
	houghinal = tf.cumsum(hough, axis = 1)

	return houghinal


def manchurian_candidates(numCand, freqIniz, image, coord):
	candidati = numpy.zeros((numCand*2,9))

	minDistance = enhancement*4

	primaFreq = freqIniz-(securbelt/2)*stepFreqRaffinato

	freqIniziale = struttura['basic_info'][0,0]['frin'][0,0][0,0]
	freqFinale = struttura['basic_info'][0,0]['frfi'][0,0][0,0]

	#QUI ANALOGO FUNZIONE CUT GD2
	#%time indexInizialewh = numpy.where(freqniu>freqIniziale)[0][0]
	#%time indexFinalewh = numpy.where(freqniu>freqFinale)[0][0]

	indexIniziale = ((freqIniziale-primaFreq)/stepFreqRaffinato).astype(numpy.int64)
	indexFinale = ((freqFinale-primaFreq)/stepFreqRaffinato+1).astype(numpy.int64)

	imageCand = image[:,indexIniziale:indexFinale].astype(numpy.int64)
	#imageCand = numpy.flip(imageCand,0)
	size = numpy.shape(imageCand)[1]
	freqniu = numpy.arange(0,size)*stepFreqRaffinato+freqIniziale

	maxPerColumn = numpy.amax(imageCand, axis = 0)
	rigaMax = numpy.argmax(imageCand, axis = 0)

	#######################

	stepFrequenzaNiu = maxPerColumn.size/numCand

	indiciFreq = numpy.arange(0,maxPerColumn.size,stepFrequenzaNiu)
	indiciFreq = numpy.append(indiciFreq, maxPerColumn.size)
	indiciFreq = numpy.round(indiciFreq).astype(numpy.int64)
	#print(indiciFreq)
	def statistics(ndArray):
		mediana = numpy.median(ndArray)
		sigmana = numpy.median(numpy.absolute(ndArray-mediana))/0.6745
		return mediana, sigmana

	stats = statistics(imageCand)
	medianaTot = stats[0]
	sigmanaTot = stats[1]
	##print(medianaTot, sigmanaTot)

	iniziali = numpy.concatenate(([indiciFreq[0]],indiciFreq[0:numCand-2],[indiciFreq[indiciFreq.size-3]]),0)

	finali = numpy.concatenate(([indiciFreq[2]],indiciFreq[3:numCand+1],[indiciFreq[indiciFreq.size-1]]),0)

	def statsPerCand(i):
		stat = statistics(maxPerColumn[iniziali[i]:finali[i]])#[0]
		return stat

	statPerCand = numpy.array(list(map(statsPerCand, numpy.arange(numCand))))
	#statPerCand = numpy.zeros((numCand,2))

	#for i in numpy.arange(numCand):
			#statPerCand[i] = statsPerCand(i)

	##print(statPerCand)
	medianaPerCand = statPerCand[:,0]
	sigmanaPerCand = statPerCand[:,1]

	#percorsoRob = ("/home/protoss/Documenti/TESI/HWITest/mioRob%d.mat" % inputFreq)
	#scipy.io.savemat(percorsoRob,{"stat": statPerCand,
										#"maxs": maxPerColumn,
										#"imax": rigaMax,
										#"imageCand": imageCand})


	filtro = numpy.where(medianaPerCand > 0)[0]
	counter = -1
	for i in filtro:
		inizio = indiciFreq[i]
		fine = indiciFreq[i+1]
		porzioneMaxPerColumn = maxPerColumn[inizio:fine]

		localMax = numpy.amax(porzioneMaxPerColumn)

		localInd = numpy.argmax(porzioneMaxPerColumn)

		if localMax > medianaPerCand[i] and localMax > medianaTot/2:
			counter = counter + 1
			index = indiciFreq[i] + localInd
			riga = rigaMax[index]

			candidati[counter,0] = freqniu[index]
			candidati[counter,1] = coord[0]
			candidati[counter,2] = coord[1]
			candidati[counter,3] = spindowns[riga]
			candidati[counter,4] = localMax
			candidati[counter,5] = (localMax-medianaPerCand[i])/sigmanaPerCand[i]
			candidati[counter,6] = coord[2]/2
			candidati[counter,7] = numpy.abs(coord[3]-coord[4])/4
			candidati[counter,8] = 1

			limite1 = numpy.amax([localInd-minDistance,0]).astype(numpy.int64)
			limite2 = numpy.amin([localInd+minDistance+1,porzioneMaxPerColumn.size]).astype(numpy.int64)
			
			porzioneMaxPerColumn[limite1:limite2] = 0
			secondLocMax = numpy.amax(porzioneMaxPerColumn)
			secondLocInd = numpy.argmax(porzioneMaxPerColumn)

			if numpy.absolute(secondLocInd-localInd) > 2 * minDistance and secondLocMax > medianaPerCand[i]:
				counter = counter + 1
				index = indiciFreq[i] + secondLocInd
				riga = rigaMax[index]

				candidati[counter,0] = freqniu[index]
				candidati[counter,1] = coord[0]
				candidati[counter,2] = coord[1]
				candidati[counter,3] = spindowns[riga]
				candidati[counter,4] = secondLocMax
				candidati[counter,5] = (secondLocMax-medianaPerCand[i])/sigmanaPerCand[i]
				candidati[counter,6] = coord[2]/2
				candidati[counter,7] = numpy.abs(coord[3]-coord[4])/4
				candidati[counter,8] = 2

	candidati[3,:]=numpy.round(candidati[3,:] / stepSpindown) * stepSpindown
	return candidati


#config = tf.ConfigProto()
#config.gpu_options.allow_growth=True
#sessione = tf.Session(config=config)

#sessione = tf.Session()

# CARICO DATI

tFft = 8192
#tFft = 4096
tObs = 9 #mesi
tObs = tObs*30*24*60*60
cands = 100
primaFreq = 108

percorsoQuad = ("quadHWI%d.mat" % primaFreq)
percorsoPatch = ("quadHWI%dEclNew.mat" % primaFreq)

#carico file per dopp corr
quadrato = scipy.io.loadmat(percorsoQuad)['quad'].astype(numpy.float64)
nPunti = quadrato.shape[0]
print(nPunti)
patch = scipy.io.loadmat(percorsoPatch)['quadratoEclNew'].astype(numpy.float64)

percorsoFile = 'dati/H/in_O2LH_01_0108_.mat'

#nPunti = 10
#nFiles = 2

allCands = numpy.zeros((1,nPunti,cands*2,9))


print(patch.shape, quadrato.shape)
print("\n\n\n")

struttura = scipy.io.loadmat(percorsoFile)['job_pack_0']
tempi = struttura['peaks'][0,0][0]#.astype(numpy.float32)
frequenze = struttura['peaks'][0,0][1]#.astype(numpy.float32)
pesi = (struttura['peaks'][0,0][4]+1)#.astype(numpy.float32)

#nb: picchi ha 0-tempi
#              1-frequenze
#              2-pesi

#headers vari
securbelt = 4000
#securbelt = 4000*3

#frequenze
stepFreq = 1/tFft
enhancement = 10
stepFreqRaffinato =  stepFreq/enhancement

#tempi
#epoca definita come mediana di tempi di tutto il run #WARNING da ridefinire con durata dati che prendo
epoca = (57722+57990)/2

#spindowns
spindownMin = -1e-9
spindownMax = 1e-10

#spindownMin = -2e-9
#spindownMax = 1e-9


stepSpindown = stepFreq/tObs

nstepSpindown = numpy.round((spindownMax-spindownMin)/stepSpindown).astype(numpy.int64)
spindowns = numpy.arange(0, nstepSpindown)
spindowns = numpy.multiply(spindowns,stepSpindown)
spindowns = numpy.add(spindowns, spindownMin)

#per doppler corr
veloc = struttura['basic_info'][0,0]['velpos'][0,0][0:3,:].astype(numpy.float64)
nTempi = struttura['basic_info'][0,0]['ntim'][0,0][0,0]
primoTempo = struttura['basic_info'][0,0]['tim0'][0,0][0,0]
indices = struttura['basic_info'][0,0]['index'][0,0][0]



#mettere un for
#start = time.time()
#for punto in numpy.arange(0,nPunti):

sessione = tf.Session()
freqInCorr,freqCorr, freqPerHough = doppcorr(punto)


nstepsFreq = numpy.ceil(securbelt+(numpy.amax(freqCorr)-numpy.amin(freqCorr) + stepFreq + 2*stepFreqRaffinato)/stepFreqRaffinato)
##print(nstepsFreq)
nColumns = tf.cast(nstepsFreq, dtype=tf.int32)

	#da qui si usa tensorflow
#definisco tutte le costanti necessarie
tempiTF = tf.constant(tempi,dtype=tf.float64)
pesiTF = tf.constant(pesi,dtype=tf.float32)
spindownsTF = tf.constant(spindowns, dtype=tf.float64)

tempiHM = tempiTF-epoca
tempiHM = ((tempiHM)*3600*24/stepFreqRaffinato)
tempiHM = tf.cast(tempiHM, tf.float64)

pesiHM = tf.reshape(pesiTF,(1,tf.size(pesiTF)))
pesiHM = pesiHM[0]

nRows = tf.constant(nstepSpindown, dtype=tf.int64)

#freqTF = noncorr()

freqTF = tf.constant(freqPerHough, dtype = tf.float64)
houghmap = inDaHough(punto,freqTF)
hough = sessione.run(houghmap)

#allCands[ifile,punto] = manchurian_candidates(cands, freqInCorr, hough, patch[punto])

#tf.reset_default_graph()
#sessione.close()

	#numCands = numpy.
#stop = time.time()


numAllCands = numpy.int64(allCands.size/9)
allCands = allCands.reshape(numAllCands,9)
nonzeri = numpy.nonzero(allCands[:,0])
finallCands = allCands[nonzeri]
finallCands = numpy.transpose(finallCands)



percorsoOut = ('candHWI%d.mat' % primaFreq)
#percorsoOut = 'provacand'


scipy.io.savemat(percorsoOut,{'cand' : finallCands})