In [2]:
import tensorflow as tf
import numpy
import scipy.io
import time

sessione = tf.Session()

# CARICO DATI

#tFft = 8192
tFft = 4096
tObs = 1/5 #mesi
tObs = tObs*30*24*60*60
nPunti = 2
cands = 100
primaFreq = 1/tFft

#percorsoDati = "dati/dati9mesi108HWI.mat"
#percorsoQuad = "quadHWI108.mat"
#percorsoPatch = "quadHWI108Ecl.mat"

#percorsoDati = "dati/dati9mesi108HWI.mat"
#percorsoQuad = ("quad%dLIL.mat" % tFft)
#percorsoPatch = ("quad%dEclNew.mat" % tFft)


percorsoDati = "/home/protoss/Documenti/TESI/wn100bkp/data/amplitude_4.mat"


#index = 9002 #signal
index = 4 #signal meno bello
#index = 5001
#index = 0

struttura1 = scipy.io.loadmat(percorsoDati)['H']
peakmap1 = struttura1[index].copy()
del struttura1

peakmapTOT = peakmap1


peaklinee = numpy.zeros(peakmapTOT.shape)

In [3]:
#spindowns
spindownMin = -10e-10
spindownMax = 9e-10
nstepSpindown = 100
stepSpindown = (spindownMax-spindownMin)/nstepSpindown
#stepSpindown = stepFreq/tObs 

#nstepSpindown = numpy.round((spindownMax-spindownMin)/stepSpindown).astype(numpy.int32)

spindowns = numpy.arange(0, nstepSpindown)
spindowns = numpy.multiply(spindowns,stepSpindown)
spindowns = numpy.add(spindowns, spindownMin)




peaklinee = numpy.zeros(peakmapTOT.shape)

peaklinee[200] = 1


colonne = numpy.arange(peaklinee.shape[1])
print(colonne.size)
righe = numpy.round(numpy.arange(105,117,(117-105)/128)).astype(int)
print(righe.size)

peaklinee[righe,colonne] = 1


peaklinee = numpy.flip(peaklinee,0)

#peaklinee[:,10] = 1


128
128

In [5]:
from matplotlib import pyplot
import matplotlib as mpl
%matplotlib notebook
mpl.rcParams['font.size'] = 21
pyplot.figure(figsize=(10, 8))
#a = pyplot.scatter(numpy.arange(tempiUnici.size),tempiUnici/30, s = 0.5)
a = pyplot.imshow(peaklinee, cmap = 'gray_r',origin = 'lower', interpolation = 'none',aspect = 0.9)
pyplot.xlabel('$x$',fontsize = 21)
pyplot.ylabel('$y$', fontsize = 21)
pyplot.title('Input space', fontsize = 24)
pyplot.savefig('simpeak.pdf', format='pdf')



In [6]:
import tensorflow as tf
import numpy
import scipy.io
import time

sessione = tf.Session()

# CARICO DATI

#tFft = 8192

percorsoDati = "/home/protoss/Documenti/TESI/wn100bkp/data/amplitude_4.mat"


#index = 9002 #signal
index = 4 #signal meno bello
#index = 5001
#index = 0

struttura1 = scipy.io.loadmat(percorsoDati)['H']
#peakmap1 = struttura1[index].copy()
#del struttura1

struttura2 = scipy.io.loadmat(percorsoDati)['L']
#peakmap2 = struttura2[index].copy()
#del struttura2

struttura3 = scipy.io.loadmat(percorsoDati)['V']
#peakmap3 = struttura3[index].copy()
#del struttura3


strutturaTOT = struttura1 + struttura2 + struttura3


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


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.float32)
	return freqNonCor, freqIniz


def inDaHough(freqHM):
	#calcola la hough per ogni step di spindown
	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)
		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.int32)
		houghInt = houghDiff[:,enhancement:nColumns]-houghDiff[:,0:nColumns - enhancement]
		houghInt = tf.concat([houghDiff[:,0:enhancement],houghInt],1)
		return houghInt

	
	hough = sliceInt()
	#hough = convInt()
	houghinal = tf.cumsum(hough, axis = 1)
	
	return houghinal
	#return houghDiff
    
    


tFft = 4096
tObs = 1/5 #mesi
tObs = tObs*30*24*60*60
#nPunti = 2
cands = 32
primaFreq = 1/tFft

#maximumCand = numpy.array()


peakmapTOT = strutturaTOT[index]
nonzeri = numpy.nonzero(peakmapTOT)
peakmapTOT[nonzeri] = 1

In [7]:
#sparsa = numpy.nonzero(peakmapTOT)
sparsa = numpy.nonzero(peaklinee)

frequenze,tempi = sparsa
tempi = tempi+1
frequenze = frequenze / tFft + 1
pesi = numpy.ones(sparsa[0].size)




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

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

#frequenze
stepFreq = 1/tFft
enhancement = 10
#enhancement = 1
stepFreqRaffinato =  stepFreq/enhancement
nstepsFreq = securbelt+(numpy.amax(frequenze)-numpy.amin(frequenze) + stepFreq + 2*stepFreqRaffinato)/stepFreqRaffinato


epoca = 55

#spindowns
spindownMin = -9e-10
spindownMax = 9e-10
nstepSpindown = 100
stepSpindown = (spindownMax-spindownMin)/nstepSpindown

spindowns = numpy.arange(0, nstepSpindown)
spindowns = numpy.multiply(spindowns,stepSpindown)
spindowns = numpy.add(spindowns, spindownMin)
print(numpy.where(spindowns == numpy.amin(numpy.absolute(spindowns))))
#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.float32) 

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

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

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

nColumns = tf.cast(nstepsFreq, dtype=tf.int32)

freq, freqIn = noncorr()

freqTF = tf.constant(freq, dtype = tf.float32)
houghmap = inDaHough(freqTF)
hough = sessione.run(houghmap)


(array([50]),)

In [10]:
from matplotlib import pyplot
import matplotlib as mpl
%matplotlib notebook
mpl.rcParams['font.size'] = 22
pyplot.figure(figsize=(12, 8))

a = pyplot.imshow(hough[:,400:-50], cmap = 'gray_r', origin = 'lower', interpolation = 'none', aspect = 30)

labelxtick = [0,50,100,150,200,250]
posxTick = [0,500,1000,1500,2000,2500]
pyplot.xticks(posxTick,labelxtick)

sdveri = (-0.1953/15)*(numpy.arange(101)-50)
print(sdveri)

#labelytick = numpy.array([-50, -25, 0, 25, 50])
sdlabel = numpy.array([0,17,33,50,66,83,100])
labelytick = sdveri[sdlabel]
labelytick = [-0.65, -0.43, -0.21, 0, 0.21, 0.43, 0.65]
posytick = sdlabel#numpy.arange(0,101,25)#[  0, 25 ,50, 75, 100]



pyplot.yticks(posytick, labelytick)

pyplot.colorbar(shrink = 0.5,aspect = 15)
pyplot.xlabel('$q$',fontsize = 22)
pyplot.ylabel('$m$',fontsize = 22)
pyplot.title('Parameters space', fontsize = 25)
pyplot.tight_layout()
pyplot.savefig('simphough.pdf', format='pdf')


[ 0.651    0.63798  0.62496  0.61194  0.59892  0.5859   0.57288  0.55986
  0.54684  0.53382  0.5208   0.50778  0.49476  0.48174  0.46872  0.4557
  0.44268  0.42966  0.41664  0.40362  0.3906   0.37758  0.36456  0.35154
  0.33852  0.3255   0.31248  0.29946  0.28644  0.27342  0.2604   0.24738
  0.23436  0.22134  0.20832  0.1953   0.18228  0.16926  0.15624  0.14322
  0.1302   0.11718  0.10416  0.09114  0.07812  0.0651   0.05208  0.03906
  0.02604  0.01302 -0.      -0.01302 -0.02604 -0.03906 -0.05208 -0.0651
 -0.07812 -0.09114 -0.10416 -0.11718 -0.1302  -0.14322 -0.15624 -0.16926
 -0.18228 -0.1953  -0.20832 -0.22134 -0.23436 -0.24738 -0.2604  -0.27342
 -0.28644 -0.29946 -0.31248 -0.3255  -0.33852 -0.35154 -0.36456 -0.37758
 -0.3906  -0.40362 -0.41664 -0.42966 -0.44268 -0.4557  -0.46872 -0.48174
 -0.49476 -0.50778 -0.5208  -0.53382 -0.54684 -0.55986 -0.57288 -0.5859
 -0.59892 -0.61194 -0.62496 -0.63798 -0.651  ]

In [8]:
colonne = numpy.arange(peakmapTOT.shape[1])
print(colonne.size)
righe = numpy.round(numpy.arange(135,151,(151-135)/128)).astype(int)
righe = numpy.flip(righe,0)
print(righe.size)

peakmapTOT[righe,colonne] = 1


128
128

In [9]:
from matplotlib import pyplot
%matplotlib notebook
pyplot.figure(figsize=(12, 8))




#peaklinee = numpy.flip(peaklinee,0)

nonzeri = numpy.nonzero(peakmapTOT)
peakmapTOT[nonzeri] = 1

a = pyplot.imshow(peakmapTOT, cmap = 'binary', origin = 'lower', interpolation = 'none', aspect = 0.9)



#sdveri = (-0.1953/15)*(numpy.arange(101)-50)
#print(sdveri)

#labelytick = numpy.array([-50, -25, 0, 25, 50])
#sdlabel = numpy.array([0,17,33,50,66,83,100])
#labelytick = sdveri[sdlabel]
#labelytick = [-0.65, -0.43, -0.21, 0, 0.21, 0.43, 0.65]

labelxtick = numpy.array([0,20.0,40.0,60.0,80.0,100.0,120.0])*8192/3600
print(labelxtick)
labelytick = numpy.array([0,50.0,100.0,150.0,200.0,250.0])
labelytick = labelytick/8192+80
print(labelytick)


pyplot.xlabel('$t$ $(hours)$')
labelxtick = [0,23,46,69,92,115,138]
posxTick = [0,20,40,60,80,100,120]
pyplot.xticks(posxTick,labelxtick)

pyplot.ylabel('$\\nu$ ($\\times 10^{-3} + 80 Hz$)' )
labelytick = [0,6,12,18,24,30]
posyTick = [0,50,100,150,200,250]
pyplot.yticks(posyTick,labelytick)
pyplot.title('Peakmap')

pyplot.show()


[  0.          45.51111111  91.02222222 136.53333333 182.04444444
 227.55555556 273.06666667]
[80.         80.00610352 80.01220703 80.01831055 80.02441406 80.03051758]

In [11]:
import tensorflow as tf
import numpy
import scipy.io
import time

sessione = tf.Session()

# CARICO DATI

#tFft = 8192

percorsoDati = "/home/protoss/Documenti/TESI/wn100bkp/data/amplitude_4.mat"


#index = 9002 #signal
index = 4 #signal meno bello
#index = 5001
#index = 0

struttura1 = scipy.io.loadmat(percorsoDati)['H']
#peakmap1 = struttura1[index].copy()
#del struttura1

struttura2 = scipy.io.loadmat(percorsoDati)['L']
#peakmap2 = struttura2[index].copy()
#del struttura2

struttura3 = scipy.io.loadmat(percorsoDati)['V']
#peakmap3 = struttura3[index].copy()
#del struttura3


strutturaTOT = struttura1 + struttura2 + struttura3


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


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.float32)
	return freqNonCor, freqIniz


def inDaHough(freqHM):
	#calcola la hough per ogni step di spindown
	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)
		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.int32)
		houghInt = houghDiff[:,enhancement:nColumns]-houghDiff[:,0:nColumns - enhancement]
		houghInt = tf.concat([houghDiff[:,0:enhancement],houghInt],1)
		return houghInt

	
	hough = sliceInt()
	#hough = convInt()
	houghinal = tf.cumsum(hough, axis = 1)
	
	return houghinal
	#return houghDiff
    
    


tFft = 4096
tObs = 1/5 #mesi
tObs = tObs*30*24*60*60
#nPunti = 2
cands = 32
primaFreq = 1/tFft

#maximumCand = numpy.array()

nIndexes = 10
#nIndexes = strutturaTOT.shape[0]
allMaxCands = numpy.zeros((nIndexes,4,5))

peakmapTOT = strutturaTOT[index]

#del struttura1
#del struttura2
#del struttura3

colonne = numpy.arange(peakmapTOT.shape[1])
print(colonne.size)
righe = numpy.round(numpy.arange(135,151,(151-135)/128)).astype(int)
righe = numpy.flip(righe,0)
print(righe.size)

peakmapTOT[righe,colonne] = 1

sparsa = numpy.nonzero(peakmapTOT)

frequenze,tempi = sparsa
tempi = tempi+1
frequenze = frequenze / tFft + 1
pesi = numpy.ones(sparsa[0].size)




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

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

#frequenze
stepFreq = 1/tFft
enhancement = 10
#enhancement = 1
stepFreqRaffinato =  stepFreq/enhancement
nstepsFreq = securbelt+(numpy.amax(frequenze)-numpy.amin(frequenze) + stepFreq + 2*stepFreqRaffinato)/stepFreqRaffinato


epoca = 55

#spindowns
spindownMin = -10e-10
spindownMax = 9e-10
nstepSpindown = 90
stepSpindown = (spindownMax-spindownMin)/nstepSpindown

spindowns = numpy.arange(0, nstepSpindown)
spindowns = numpy.multiply(spindowns,stepSpindown)
spindowns = numpy.add(spindowns, spindownMin)

#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.float32) 

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

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

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

nColumns = tf.cast(nstepsFreq, dtype=tf.int32)

freq, freqIn = noncorr()

freqTF = tf.constant(freq, dtype = tf.float32)
houghmap = inDaHough(freqTF)
hough = sessione.run(houghmap)


from matplotlib import pyplot
%matplotlib notebook
pyplot.figure(figsize=(12, 8))



a = pyplot.imshow(hough[:,1000:-1000], cmap = 'gray_r', origin = 'lower', interpolation = 'none', aspect = 30)


labelxtick = [0,6,12,18,24,30]
posxTick = [0,500,1000,1500,2000,2500]
#pyplot.xticks(posxTick,labelxtick)

labelytick = numpy.arange(-8,8.1,2)
print(labelytick)
posyTick = [7,17,27,37,47, 57,67,77,87]

#pyplot.yticks(posytick,labelytick)
pyplot.xlabel('$\\nu_0$ ($\\times 10^{-3} + 80 Hz$)')
pyplot.ylabel('$\dot{\\nu} \;(10^{-10}\:Hz/s)$')

pyplot.xticks(posxTick,labelxtick)
pyplot.yticks(posyTick,labelytick)

pyplot.colorbar(shrink = 0.5,aspect = 15)
pyplot.title('Frequency-Hough map')
pyplot.show()


128
128
[-8. -6. -4. -2.  0.  2.  4.  6.  8.]

In [35]:
zeri = numpy.where(hough[:,1000:-1000] == 0)
zeri


Out[35]:
(array([86, 86]), array([2529, 2530]))

In [59]:
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D #<-- Note the capitalization! 

%matplotlib notebook

indici = numpy.nonzero(hough[:,1000:-1000])
righe = indici[0]
colonne = indici[1]
righe = numpy.concatenate([righe,[86,86]])
colonne = numpy.concatenate([colonne,[2529,2530]])
valori = numpy.ravel(hough[:,1000:-1000])

print(righe.size, colonne.size, valori.size)



fig = plt.figure()

ax = Axes3D(fig) #<-- Note the difference from your original code...

# Plot the values
ax.scatter(righe, colonne, valori, c = valori, marker='o', alpha = 1)


228780 228780 228780
Out[59]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fad22bfb4e0>

In [26]:
from scipy import sparse

sparsh = sparse.csr_matrix(hough[:,1000:-1000])

In [30]:
print(sparsh.shape)


(90, 2542)

In [ ]: