unimatr = numpy.ones((10,10))
#unimatr
duimatr = unimatr*2
#duimatr

uniarray = numpy.ones((10,1))
#uniarray

triarray = uniarray*3

scalarray = numpy.arange(10)
scalarray = scalarray.reshape(10,1)

#NB fare il reshape da orizzontale a verticale è come se aggiungesse
#una dimensione all'array facendolo diventare un ndarray
#(prima era un array semplice, poi diventa un array (x,1), quindi puoi fare trasposto)
#NB NUMPY NON FA TRASPOSTO DI ARRAY SEMPLICE!
#scalarray
scalarray.T

ramatricia = numpy.random.randint(2, size=36).reshape((6,6))
ramatricia2 = numpy.random.randint(2, size=36).reshape((6,6))

#WARNING questa operazione moltiplica elemento per elemento
#se l'oggetto è di dimensione inferiore moltiplica ogni riga/colonna
# o matrice verticale/orizzontale a seconda della forma dell'oggetto

duimatr*scalarray
#duimatr*scalarray.T
#duimatr*duimatr
ramatricia*ramatricia2

#numpy dot invece fa prodotto matriciale righe per colonne

numpy.dot(duimatr,scalarray)
#numpy.dot(duimatr,duimatr)
numpy.dot(ramatricia2,ramatricia)

duimatr + scalarray

"""
Explanation: Prove manipolazioni array
End of explanation
"""


scalarray = numpy.arange(10)
uniarray = numpy.ones(10)

matricia = numpy.outer(scalarray, uniarray)
matricia

tensorio = numpy.outer(matricia,scalarray).reshape(10,10,10)
tensorio
# metodo di creazione array nd (numpy.ndarray)

"""
Explanation: Prove creazione matrice 3D con prodotti esterni
End of explanation
"""


tensorio = numpy.ones(1000).reshape(10,10,10)
tensorio
# metodo di creazione array nd (numpy.ndarray)
#altro metodo è con comando diretto
#tensorio = numpy.ndarray((3,3,3), dtype = int, buffer=numpy.arange(30))
#potrebbe essere utile con la matrice sparsa della peakmap, anche se difficilmente è maneggiabile come matrice densa
#oppure

# HO FINALMENTE SCOPERTO COME SI METTE IL DTYPE COME SI DEVE!! con "numpy.float32"!
#tensorio = numpy.zeros((3,3,3), dtype = numpy.float32)
#tensorio.dtype
#tensorio


scalarray = numpy.arange(10)
uniarray = numpy.ones(10)
scalamatricia = numpy.outer(scalarray,scalarray)
#scalamatricia


tensorio * 2
tensorio + 2
tensorio + scalamatricia
%time tensorio + scalarray
%time tensorio.__add__(scalarray)
#danno stesso risultato con tempi paragonabili


"""
Explanation: Prove manipolazione matrici 3D numpy
End of explanation
"""


from scipy import sparse


ramatricia = numpy.random.randint(2, size=25).reshape((5,5))
ramatricia

#efficiente per colonne
#sparsamatricia = sparse.csc_matrix(ramatricia)
#print(sparsamatricia)

#per righe
sparsamatricia = sparse.csr_matrix(ramatricia)
print(sparsamatricia)

sparsamatricia.toarray()

righe = numpy.array([0,0,0,1,2,3,3,4])
colonne = numpy.array([0,0,4,2,1,4,3,0])
valori = numpy.ones(righe.size)
sparsamatricia = sparse.coo_matrix((valori, (righe,colonne)))

print(sparsamatricia)

sparsamatricia.toarray()

"""
Explanation: Prove matrici sparse
End of explanation
"""


#vari modi per fare prodotti di matrici (con somma con operatore + è lo stesso)
densamatricia = sparsamatricia.toarray()

#densa-densa
prodottoPerElementiDD = densamatricia*densamatricia
prodottoMatricialeDD = numpy.dot(densamatricia, densamatricia)

#sparsa-densa
prodottoMatricialeSD = sparsamatricia*densamatricia
prodottoMatricialeSD2 = sparsamatricia.dot(densamatricia)

#sparsa-sparsa
prodottoMatricialeSS = sparsamatricia*sparsamatricia
prodottoMatricialeSS2 = sparsamatricia.dot(sparsamatricia)

# "SPARSA".dot("SPARSA O DENSA") FA PRODOTTO MATRICIALE
# "SPARSA * SPARSA" FA PRODOTTO MATRICIALE


prodottoMatricialeDD - prodottoMatricialeSS
#nb somme e sottrazioni tra matrici sparse e dense sono ok
# prodotto matriciale tra densa e sparsa funziona come sparsa e sparsa

"""
Explanation: Prodotto di matrici
Prodotti interni
Considera di avere 2 matrici, a e b, in forma numpy array:

a*b fa il prodotto elemento per elemento (solo se a e b hanno stessa dimensione)
numpy.dot(a,b) fa il prodotto matriciale righe per colonne

Ora considera di avere 2 matrici, a e b, in forma di scipy.sparse:

a*b fa il prodotto matriciale righe per colonne
numpy.dot(a,b) non funziona per nulla
a.dot(b) fa il prodotto matriciale righe per colonne
End of explanation
"""


densarray = numpy.array(["a","b"],dtype = object)
densarray2 = numpy.array(["c","d"],dtype = object)

numpy.outer(densarray,[1,2])

densamatricia = numpy.array([[1,2],[3,4]])
densamatricia2 = numpy.array([["a","b"],["c","d"]], dtype = object)
numpy.outer(densamatricia2,densamatricia).reshape(4,2,2)

densarray1 = numpy.array([0,2])
densarray2 = numpy.array([5,0])
densamatricia = numpy.array([[1,2],[3,4]])
densamatricia2 = numpy.array([[0,2],[5,0]])

nrighe = 2
ncolonne = 2
npiani = 4
prodottoEstDD = numpy.outer(densamatricia,densamatricia2).reshape(npiani,ncolonne,nrighe)
#prodottoEstDD
#prodottoEstDD = numpy.dstack((prodottoEstDD[0,:],prodottoEstDD[1,:]))

prodottoEstDD


sparsarray1 = sparse.csr_matrix(densarray1)
sparsarray2 = sparse.csr_matrix(densarray2)
sparsamatricia = sparse.csr_matrix(densamatricia)
sparsamatricia2 = sparse.csr_matrix(densamatricia2)

prodottoEstSS = sparse.kron(sparsamatricia,sparsamatricia2).toarray()

prodottoEstSD = sparse.kron(sparsamatricia,densamatricia2).toarray()
prodottoEstSD




#prove prodotti esterni
# numpy.outer
# scipy.sparse.kron

#densa-densa
prodottoEsternoDD = numpy.outer(densamatricia,densamatricia)

#sparsa-densa
prodottoEsternoSD = sparse.kron(sparsamatricia,densamatricia)

#sparsa-sparsa
prodottoEsternoSS = sparse.kron(sparsamatricia,sparsamatricia)


prodottoEsternoDD-prodottoEsternoSS

# altre prove di prodotti esterni
rarray1 = numpy.random.randint(2, size=4)
rarray2 = numpy.random.randint(2, size=4)
print(rarray1,rarray2)
ramatricia = numpy.outer(rarray1,rarray2)
unimatricia = numpy.ones((4,4)).astype(int)
#ramatricia2 = rarray1 * rarray2.T
print(ramatricia,unimatricia)
#print(ramatricia)
#print("eppoi")
#print(ramatricia2)

#sparsarray = sparse.csr_matrix(rarray1)
#print(sparsarray)

#ramatricia2 = 

#il mio caso problematico è che ho una matrice di cui so tutti gli elementi non zero,
#so quante righe ho (i tempi), ma non so quante colonne di freq ho
randomcolonne = numpy.random.randint(10)+1
ramatricia = numpy.random.randint(2, size=10*randomcolonne).reshape((10,randomcolonne))
print(ramatricia.shape)
#ramatricia
nonzeri = numpy.nonzero(ramatricia)
ndati = len(nonzeri[0])
ndati
ramatricia

#ora cerco di fare la matrice sparsa
print(ndati)
dati = numpy.ones(2*ndati).reshape(ndati,2)
dati
coordinateRighe = nonzeri[0]
coordinateColonne = nonzeri[1]
sparsamatricia = sparse.coo_matrix((dati,(coordinateRighe,coordinateColonne)))
densamatricia = sparsamatricia.toarray()
densamatricia

"""
Explanation: Prodotti esterni
End of explanation
"""


matrice = numpy.arange(30).reshape(10,3)
matrice

righe = numpy.array([1,0,1,1])
colonne = numpy.array([2,0,2,2])
pesi = numpy.array([100,200,300,10])
print(righe,colonne)

matrice[righe,colonne]


matrice[righe,colonne] = (matrice[righe,colonne] + numpy.array([100,200,300,10]))
matrice

%matplotlib inline
a = pyplot.imshow(matrice)

numpy.add.at(matrice, [righe,colonne],pesi)
matrice

%matplotlib inline
a = pyplot.imshow(matrice)

matr

"""
Explanation: Provo a passare operazioni a array con array di coordinate
End of explanation
"""


from matplotlib import pyplot
%matplotlib inline


##AL MOMENTO INUTILE, NON COMPILARE
x = numpy.random.randint(10,size = 10)
y = numpy.random.randint(10,size = 10)
pyplot.scatter(x,y, s = 5)
#nb imshow si può fare solo con un 2d array

#visualizzazione di una matrice, solo matrici dense a quanto pare
a = pyplot.imshow(densamatricia)
#a = pyplot.imshow(sparsamatricia)
#c = pyplot.matshow(densamatricia)


#spy invece funziona anche per le sparse!
pyplot.spy(sparsamatricia,precision=0.01, marker = ".", markersize=10)

#in alternativa, scatterplot delle coordinate dal dataframe
b = pyplot.scatter(coordinateColonne,coordinateRighe, s = 2)

import seaborn
%matplotlib inline


sbRegplot = seaborn.regplot(x=coordinateRighe, y=coordinateColonne, color="g", fit_reg=False)

import pandas

coordinateRighe = coordinateRighe.reshape(len(coordinateRighe),1)
coordinateColonne = coordinateColonne.reshape(len(coordinateColonne),1)
#print([coordinateRighe,coordinateColonne])
coordinate = numpy.concatenate((coordinateRighe,coordinateColonne),axis = 1)
coordinate


tabella = pandas.DataFrame(coordinate)
tabella.columns = ["righe", "colonne"]


sbPlmplot = seaborn.lmplot(x = "righe", y = "colonne", data = tabella, fit_reg=False)



"""
Explanation: Prove plots
End of explanation
"""


import numpy
from scipy import sparse
import multiprocessing
from matplotlib import pyplot

#first i build a matrix of some x positions vs time datas in a sparse format
matrix = numpy.random.randint(2, size = 100).astype(float).reshape(10,10)
x = numpy.nonzero(matrix)[0]
times = numpy.nonzero(matrix)[1]
weights = numpy.random.rand(x.size)



import scipy.io

mint = numpy.amin(times)
maxt = numpy.amax(times)

scipy.io.savemat('debugExamples/numpy.mat',{
    'matrix':matrix, 
    'x':x, 
    'times':times, 
    'weights':weights,
    'mint':mint,
    'maxt':maxt,
    
})

times

#then i define an array of y positions
nStepsY = 5
y = numpy.arange(1,nStepsY+1)

# provo a iterare
# VERSIONE CON HACK CON SPARSE verificato viene uguale a tutti gli altri metodi più semplici che ho provato
# ma ha problemi con parallelizzazione

nRows = nStepsY
nColumns = 80
y = numpy.arange(1,nStepsY+1)
image = numpy.zeros((nRows, nColumns))
def itermatrix(ithStep):
    yTimed = y[ithStep]*times
    positions = (numpy.round(x-yTimed)+50).astype(int)

    fakeRow = numpy.zeros(positions.size)
    matrix = sparse.coo_matrix((weights, (fakeRow, positions))).todense()
    matrix = numpy.ravel(matrix)
    missColumns = (nColumns-matrix.size)
    zeros = numpy.zeros(missColumns)
    matrix = numpy.concatenate((matrix, zeros))
    return matrix

#for i in numpy.arange(nStepsY):
#    image[i] = itermatrix(i)

#or
imageSparsed = list(map(itermatrix, range(nStepsY)))
imageSparsed = numpy.array(imageSparsed)
scipy.io.savemat('debugExamples/numpyResult.mat', {'imageSparsed':imageSparsed})    
a = pyplot.imshow(imageSparsed, aspect = 10)
pyplot.show()

import numpy
from scipy import sparse
import multiprocessing
from matplotlib import pyplot

#first i build a matrix of some x positions vs time datas in a sparse format
matrix = numpy.random.randint(2, size = 100).astype(float).reshape(10,10)
times = numpy.nonzero(matrix)[0]
freqs = numpy.nonzero(matrix)[1]
weights = numpy.random.rand(times.size)

#then i define an array of y positions
nStepsSpindowns = 5
spindowns = numpy.arange(1,nStepsSpindowns+1)


#PROVA CON BINCOUNT

def mapIt(ithStep):
    ncolumns = 80
    image = numpy.zeros(ncolumns)

    sdTimed = spindowns[ithStep]*times
    positions = (numpy.round(freqs-sdTimed)+50).astype(int)

    values = numpy.bincount(positions,weights)
    values = values[numpy.nonzero(values)]
    positions = numpy.unique(positions)
    image[positions] = values
    return image


%time    imageMapped = list(map(mapIt, range(nStepsSpindowns)))
imageMapped = numpy.array(imageMapped)

%matplotlib inline
a = pyplot.imshow(imageMapped, aspect = 10)

# qui provo fully vectorial
def fullmatrix(nRows, nColumns):
    spindowns = numpy.arange(1,nStepsSpindowns+1)
    image = numpy.zeros((nRows, nColumns))

    sdTimed = numpy.outer(spindowns,times)
    freqs3d = numpy.outer(numpy.ones(nStepsSpindowns),freqs)
    weights3d = numpy.outer(numpy.ones(nStepsSpindowns),weights)
    spindowns3d = numpy.outer(spindowns,numpy.ones(times.size))
    positions = (numpy.round(freqs3d-sdTimed)+50).astype(int)

    matrix = sparse.coo_matrix((numpy.ravel(weights3d), (numpy.ravel(spindowns3d), numpy.ravel(positions)))).todense()
    return matrix

%time image = fullmatrix(nStepsSpindowns, 80)
a = pyplot.imshow(image, aspect = 10)
pyplot.show()

"""
Explanation: Un esempio semplice del mio problema
End of explanation
"""


#confronto con codice ORIGINALE in matlab
immagineOrig = scipy.io.loadmat('debugExamples/dbOrigResult.mat')['binh_df0']
a = pyplot.imshow(immagineOrig[:,0:80], aspect = 10)
pyplot.show()

#PROVA CON BINCOUNT

def mapIt(ithStep):
    ncolumns = 80
    image = numpy.zeros(ncolumns)

    yTimed = y[ithStep]*times
    positions = (numpy.round(x-yTimed)+50).astype(int)

    values = numpy.bincount(positions,weights)
    values = values[numpy.nonzero(values)]
    positions = numpy.unique(positions)
    image[positions] = values
    return image


%time    imageMapped = list(map(mapIt, range(nStepsY)))
imageMapped = numpy.array(imageMapped)

%matplotlib inline
a = pyplot.imshow(imageMapped, aspect = 10)

# qui provo con vettorializzazione di numpy (apply along axis)
nrows = nStepsY
ncolumns = 80
matrix = numpy.zeros(nrows*ncolumns).reshape(nrows,ncolumns)

def applyIt(image):
    ithStep = 1
    image = numpy.zeros(ncolumns)

    yTimed = y[ithStep]*times
    positions = (numpy.round(x-yTimed)+50).astype(int)
    #print(positions)
    values = numpy.bincount(positions,weights)
    values = values[numpy.nonzero(values)]
    positions = numpy.unique(positions)
    image[positions] = values
    
    return image


imageApplied = numpy.apply_along_axis(applyIt,1,matrix)
a = pyplot.imshow(imageApplied, aspect = 10)

# qui provo fully vectorial
def fullmatrix(nRows, nColumns):
    y = numpy.arange(1,nStepsY+1)
    image = numpy.zeros((nRows, nColumns))

    yTimed = numpy.outer(y,times)
    x3d = numpy.outer(numpy.ones(nStepsY),x)
    weights3d = numpy.outer(numpy.ones(nStepsY),weights)
    y3d = numpy.outer(y,numpy.ones(x.size))
    positions = (numpy.round(x3d-yTimed)+50).astype(int)

    matrix = sparse.coo_matrix((numpy.ravel(weights3d), (numpy.ravel(y3d), numpy.ravel(positions)))).todense()
    return matrix

%time image = fullmatrix(nStepsY, 80)
a = pyplot.imshow(image, aspect = 10)
pyplot.show()

imageMapped = list(map(itermatrix, range(nStepsY)))
imageMapped = numpy.array(imageMapped)
a = pyplot.imshow(imageMapped, aspect = 10)
pyplot.show()

# prova con numpy.put

nStepsY = 5

def mapIt(ithStep):
    ncolumns = 80
    image = numpy.zeros(ncolumns)

    yTimed = y[ithStep]*times
    positions = (numpy.round(x-yTimed)+50).astype(int)

    values = numpy.bincount(positions,weights)
    values = values[numpy.nonzero(values)]
    positions = numpy.unique(positions)
    image[positions] = values
    return image


%time    imagePutted = list(map(mapIt, range(nStepsY)))
imagePutted = numpy.array(imagePutted)

%matplotlib inline
a = pyplot.imshow(image, aspect = 10)
pyplot.show()

"""
Explanation: Confronti Debug!
End of explanation
"""


