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

sessione = tf.Session()

# CARICO DATI

tFft = 8192
inputFreq = 187

if inputFreq > 128:
	tFft = 4096
print(tFft)
#tFft = 4096
tObs = 9 #mesi
tObs = tObs*30*24*60*60
nPunti = 2
cands = 100



percorsoDati = ("dati/dati9mesi%dHWI.mat"% inputFreq)
percorsoQuad = ("quadHWI%d.mat"% inputFreq)
percorsoPatch = ("quadHWI%dEcl.mat"% inputFreq)
percorsoOut = ("/home/protoss/Documenti/TESI/HWITest/mio%d.mat" % inputFreq)

#percorsoQuad = ("quad%dLIL.mat" % tFft)
#percorsoPatch = ("quad%dEclNew.mat" % tFft)



#carico file dati
quadrato = scipy.io.loadmat(percorsoQuad)['quad'].astype(numpy.float64)
patch = scipy.io.loadmat(percorsoPatch)['quadratoEclNew'].astype(numpy.float64)

print(patch)

struttura = scipy.io.loadmat(percorsoDati)['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)

#print(tempi.size,frequenze.size)

#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

if inputFreq > 128:
	securbelt = 4000*2
	spindownMin = -9e-9
	spindownMax = -8e-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)

#LO STEP DI SPINDOWN LO SI DEFINISCE TRAMITE LA BANDA IN FREQUENZA DELLA PEAKMAP (in hz) 
#SU TEMPO DI OSSERVAZIONE (in sec):
#STEPFREQ/TOBS! DA QUI SCEGLIAMO QUALE SD MASSIMO E MINIMO TENERE

# COME VALORE MASSIMO TENERE SUI 1*-10^-9
# COME MASSIMO ANCHE CIRCA 1*+10^-10
# CONSIDERARE EVENTUALE SORGENTE NELLA SCELTA DI INTERVALLO SPINDOWN

#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]

#SPLITTO DOPP CORR E HOUGH TRANS

# 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)
		
		#valori = tf.bincount(appoggio,weights=pesiTF)

		#zeriDopo = tf.zeros([nColumns - tf.size(valori)], dtype=tf.float64)
		
		#riga = tf.concat([valori,zeriDopo],0, name = "rigadihough")
		#return riga


	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)
		
		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):
	minDistance = enhancement*4

	candidati = numpy.zeros((9,numCand*2))

	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]
	start = time.time()

	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 = 0
	for i in filtro:
		inizio = indiciFreq[i]
		fine = indiciFreq[i+1]
		porzioneMaxPerColumn = maxPerColumn[inizio:fine]
		
		localMax = numpy.amax(porzioneMaxPerColumn)
		
		localInd = numpy.argmax(porzioneMaxPerColumn)
		
		if i == 1:
			print(inizio, fine)
			print(porzioneMaxPerColumn, porzioneMaxPerColumn.size)
			print(localMax, localInd)
			print(medianaPerCand[i])
			print(medianaTot/2)
		
		
		
		
		if localMax > medianaPerCand[i] and localMax > medianaTot/2:
			counter = counter + 1
			index = indiciFreq[i] + localInd
			riga = rigaMax[index]
			
			candidati[0,counter] = freqniu[index]
			candidati[1,counter] = coord[0]
			candidati[2,counter] = coord[1]
			candidati[3,counter] = spindowns[riga]
			candidati[4,counter] = localMax
			candidati[5,counter] = (localMax-medianaPerCand[i])/sigmanaPerCand[i]
			candidati[6,counter] = coord[2]/2
			candidati[7,counter] = numpy.abs(coord[3]-coord[4])/4
			candidati[8,counter] = 1

			limite1 = numpy.amax([localInd-minDistance,1]).astype(numpy.int64)
			limite2 = numpy.amin([localInd+minDistance,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[0,counter] = freqniu[index]
				candidati[1,counter] = coord[0]
				candidati[2,counter] = coord[1]
				candidati[3,counter] = spindowns[riga]
				candidati[4,counter] = secondLocMax
				candidati[5,counter] = (secondLocMax-medianaPerCand[i])/sigmanaPerCand[i]
				candidati[6,counter] = coord[2]/2
				candidati[7,counter] = numpy.abs(coord[3]-coord[4])/4
				candidati[8,counter] = 2

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


    
    
#da qui si usa tensorflow
#definisco tutte le costanti necessarie
tempiTF = tf.constant(tempi,dtype=tf.float64)
pesiTF = tf.constant(pesi,dtype=tf.float64)
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)
#problema! num step freq cambia a seconda della correzione doppler
#perché freq min e freq max possono variare e lo step lo si lascia uguale
#posso andarci in 2 modi: uno è tagliando a 96000 tutto
#uno è mettendo un po' di zeri prima e dopo, cercando con la doppler corr quale è la max assoluta
#e quale è la min assoluta


#freqTF = noncorr()

#mettere un for
start = time.time()
for punto in numpy.arange(0,nPunti-1):
	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)
	
	freqTF = tf.constant(freqPerHough, dtype = tf.float64)
	houghmap = inDaHough(punto,freqTF)
	hough = sessione.run(houghmap)
	candidati = manchurian_candidates(cands, freqInCorr, hough, patch[punto])
	nonzeri = numpy.nonzero(candidati[0])
	finalCand = candidati[:,nonzeri]
	
	
stop = time.time()
print(stop-start)
#print(freqCorr)

scipy.io.savemat(percorsoOut,{"freqCorr": freqCorr,
						"hough": hough,
						"candidati": finalCand
						})


from matplotlib import pyplot


pyplot.figure(figsize=(10, 8))
##posxTick = numpy.arange(5551, 89454, round((89454-5551)/10))
#posxTick = numpy.arange(2330, 210300, round((210300-2330)/50))
##labelxTick = numpy.arange(0,1.1,0.1)+inputFreq
#labelxTick = numpy.arange(0,5.1,0.1)+inputFreq
#pyplot.xticks(posxTick,labelxTick)


#posyTick = numpy.arange(11)*nstepSpindown/10
#labelyTick = numpy.arange(spindownMin, spindownMax, stepSpindown*nstepSpindown/10)
#pyplot.yticks(posyTick,labelyTick)





a = pyplot.imshow(hough, aspect = 400)
pyplot.colorbar(shrink = 1,aspect = 10)

pyplot.show()