TUTORIAL PROFILE PLOT

In questo tutorial viene mostrato come costruire in modo rapido ed indolore un profile plot usando il pitone.

Il punto di partenza è l'esempio di Valerio da cui prenderò anche ispirazione per i dati da usare.

import numpy as np
from matplotlib import pyplot as plt

Genero e plotto dei dati casuali. Utilizzo un differente binnaggio per poter seguire al meglio quello che succede, e non trovarmi a dover riflettere se lavorare sulla matrice hh[0] o sulla sua trasposta. Lungo la prima dimensione (ovvero muovendosi lungo le righe) viene riportata la suddivisione in x, che nel plot viene mostrata in orizzontale

N = 100000
xVect = np.random.rand(N)*2*np.pi # flat distribution between 0 and 2pi
yVect = np.sin(xVect) + np.random.normal(0, 0.1, N) +  np.random.normal(0, 0.3, N)*xVect/5 # sin wave + some spread which increases with x 

# Numero di bin
Nx = 50
Ny = 60

fig, ax = plt.subplots()
fig.set_size_inches(11, 6)

mat, binsx, binsy, hh = ax.hist2d(xVect, yVect, bins = (Nx, Ny), cmap='coolwarm')

fig.colorbar(hh, ax = ax)
plt.show()

Definisco una comoda funzione che crea il profile plot.

Ricordiamo che i for in python sono il male assoluto quando si itera sull'intero set di dati. Iterare sui bin dell'istogramma sfruttando la vettorializzazione, è in generale una cosa lecita e comoda (come testimoniano i tempi di esecuzione)

def myProfile(xVect, yVect, nbin):
    h, bins = np.histogram(xVect, bins = nbin)
    binc = bins[:-1] + (bins[1]- bins[0])/2
    
    out_vect = np.array([np.mean(yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbin)])
    out_err =  np.array([np.std (yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbin)])
    
    return binc, out_vect, out_err

%%time
fig, ax = plt.subplots()
fig.set_size_inches(11, 6)

mat, binsx, binsy, hh = ax.hist2d(xVect, yVect, bins = (Nx, Ny), cmap='coolwarm')


%time binc, out_vect, out_err = myProfile(xVect, yVect, Nx)
ax.errorbar(binc, out_vect, yerr = out_err, ls = ":", c = "k", marker = ".")

fig.colorbar(hh, ax = ax)
plt.show()

CPU times: total: 62.5 ms
Wall time: 101 ms

CPU times: total: 422 ms
Wall time: 620 ms

Cosa succede se l'istogramma 2D ha un range verticale diverso da quello dei dati?

Infatti, la funzione sopra definisce per ciascuna fetta verticale la media sui dati, senza avere memoria del range verticale dell'istogramma

%%time
fig, ax = plt.subplots()
fig.set_size_inches(11, 6)

mat, binsx, binsy, hh = ax.hist2d(xVect, yVect, bins = (Nx, Ny), cmap='coolwarm', range = ((0, 2*np.pi),(-.5, .5)))


%time binc, out_vect, out_err = myProfile(xVect, yVect, Nx)
ax.errorbar(binc, out_vect, yerr = out_err, ls = ":", c = "k", marker = ".")

fig.colorbar(hh, ax = ax)
plt.show()

CPU times: total: 15.6 ms
Wall time: 15.4 ms

CPU times: total: 250 ms
Wall time: 290 ms

Definiamo quindi una funzione che tenga traccia anche del range verticale

def myProfile(xVect, yVect, nbinx, rangex, rangey):
    h, bins = np.histogram(xVect, bins = nbinx, range = rangex)
    binc = bins[:-1] + (bins[1]- bins[0])/2
    
    log = (yVect >= rangey[0]) & (yVect <= rangey[1])
    xVect = xVect[log]
    yVect = yVect[log]
    
    out_vect = np.array([np.nanmean(yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbinx)])
    out_err =  np.array([np.nanstd (yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbinx)])
    
    return binc, out_vect, out_err

Indipendentemente dal significato che possa avere, questa funzione fa il profile plot SOLO DI QUELLO CHE SI VEDE NELLA FIGURA, che è il comportamento standard di paw

%%time
rangex = (0, 2*np.pi)
rangey = (-.5, .5)

fig, ax = plt.subplots()
fig.set_size_inches(11, 6)

mat, binsx, binsy, hh = ax.hist2d(xVect, yVect, bins = (Nx, Ny), cmap='coolwarm', range = (rangex, rangey))


%time binc, out_vect, out_err = myProfile(xVect, yVect, Nx, rangex, rangey)
ax.errorbar(binc, out_vect, yerr = out_err, ls = ":", c = "k", marker = ".")

fig.colorbar(hh, ax = ax)
plt.show()

C:\Users\steca\AppData\Local\Temp\ipykernel_21000\3407349314.py:9: RuntimeWarning: Mean of empty slice
  out_vect = np.array([np.nanmean(yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbinx)])
c:\Users\steca\anaconda3\envs\AnaTesiLM\Lib\site-packages\numpy\lib\nanfunctions.py:1879: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,

CPU times: total: 31.2 ms
Wall time: 22.9 ms

CPU times: total: 281 ms
Wall time: 291 ms

Definiamo infine una funzione che faccia in un colpo solo istogramma 2D e profile plot

def Histo2DProfile(xVect, yVect, nbinx, nbiny, rangex, rangey, **kwargs):
    
    fig, ax = plt.subplots()
    fig.set_size_inches(11, 6)

    # Hist 2d
    mat, binsx, binsy, hh = ax.hist2d(xVect, yVect, bins = (nbinx, nbiny), range = (rangex, rangey), **kwargs)
    
    
    # Plug here the old function
    h, bins = np.histogram(xVect, bins = nbinx, range = rangex)
    binc = bins[:-1] + (bins[1]- bins[0])/2
    
    log = (yVect >= rangey[0]) & (yVect <= rangey[1])
    xVect = xVect[log]
    yVect = yVect[log]
    
    out_vect = np.array([np.nanmean(yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbinx)])
    out_err =  np.array([np.nanstd (yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbinx)])
    
    # Plot profile plot
    ax.errorbar(binc, out_vect, yerr = out_err, ls = ":", c = "k", marker = ".")

    fig.colorbar(hh, ax = ax)
    plt.show()

# Proviamo la funzione
Histo2DProfile(xVect, yVect, Nx, Ny, rangex, rangey, cmap='coolwarm')

C:\Users\steca\AppData\Local\Temp\ipykernel_21000\614021229.py:17: RuntimeWarning: Mean of empty slice
  out_vect = np.array([np.nanmean(yVect[(xVect >= bins[i]) & (xVect < bins[i+1])]) for i in range(nbinx)])