kNNpy.Auxiliary.PeakStatistics

  1import numpy as np
  2import matplotlib.pyplot as plt
  3from pylab import *
  4import scipy
  5import healpy as hp
  6
  7def PeakCurves(DensityFields=[], Nreals=10, MaxThreshold=16, Nthresh=101, Type=0, Plot=1, LogScale=1, CosmoLabels=['null']):
  8    '''
  9    Gives the peak curves for the given cosmologies' square (over)density fields.
 10    
 11    Parameters
 12    ----------
 13    DensityFields: float array of shape ''(nCosmo, NofRealizations, XdfDim, YdfDim)''
 14        4D Array of the 2D (over)density fields of the various cosmologies to be compared. The array should of shape (nCosmo, NofRealizations, XdfDim, YdfDim) 
 15        where 'nCosmo' is the number of cosmologies to be compared, 'NofRealizations' is the number of realizations of each input cosmology 
 16        (NofRealizations >= Nreals), and 'XdfDim' and 'YdfDim' are the dimensions of the 2D density fields in pixels.
 17        Example: np.array(DensityFields).shape = (3, 10, 512, 512) - 3 cosmologies containing 10 realisations each of (512x512) pixel 2D density fields.
 18    Nreals: int
 19        Number of realisations of the density fields to be used for the peak curves' statistics. Naturally, cannot be larger than the inherent number 
 20        of realisations of the 2D density fields contained within the (.npy) data files input, i.e. in the above example, Nreals <=10. Any non-int values 
 21        are type-cast to int.
 22    MaxThreshold: float
 23        Maximum overdensity threshold for which the peak values are to be plotted.
 24    Nthresh: int
 25        Number of overdensity threshold (X-axis) values to be computed in the closed interval [-1,MaxThreshold]. Any non-int values are type-cast to int.
 26    Type: int: 0 or 1, optional
 27        Type of peak curve to be plotted - 0 for raw peak curve plot; 1 for peak curves normalized by the first input cosmology's peak curve.
 28    Plot: int: 0 or 1, optional
 29        1 to plot the output and return results, 0 to skip plotting and only return results. Any values other than 1 will skip plotting.
 30    LogScale: int: 0 or 1, optional
 31        1 to plot the peak curves on log Y-axis, 0 to not use log Y-axis. Any values other than 1 will not output the log scale.
 32    CosmoLabels: str array of shape ''(nCosmo)''
 33        List of the names/labels to be assigned to the respective input cosmologies. Must have length equal to the number of cosmologies input ('nCosmo').
 34        Example: ['CDM', 'SIDM', 'FDM']
 35            
 36    Returns
 37    -------
 38    thresh, tmean, tstddev: 3 numpy arrays containing the threshold values, mean peak values and std. dev. of the peak values
 39        For (Type = 0): Peak curve plot of the various input comologies' density fields.
 40        For (Type = 1): Peak curve plot of the various input comologies' density fields normalized by the first input cosmology's density field.
 41        In both cases, the threshold values array (X-axis, 1D) 'thresh', the peaks array (Y-axis, 2D) 'tmean' containing the number of peaks corresponding to 
 42        the thresholds array for each input cosmology and their corresponding standard deviations (error, 2D) 'tstddev' are also returned.
 43
 44    Raises
 45    ------
 46    ValueError
 47        If 'Type' is not 0 or 1.
 48    ValueError
 49        If 'Nreals' is lesser than 1.
 50    ValueError
 51        If 'MaxThreshold' is lesser than or equal to (-1).
 52    ValueError
 53        If the plot needs to be output and the number of labels ('CosmoLabels') is not equal to the number of input cosmologies.
 54    '''
 55    if((Type!=0) and (Type!=1)):
 56        print("ERROR: Incorrect output type requested. Please use either 0 or 1 for variable 'Type'.")
 57    elif(Nreals<=0):
 58        print("ERROR: Invalid number of realizations entered. Please enter a positive integer value.")
 59    elif(MaxThreshold<=-1):
 60        print("ERROR: Please enter a maximum overdensity threshold value >(-1).")
 61    elif((Plot==1) and (len(CosmoLabels)!=len(DensityFields))):
 62        print("ERROR: Please enter the labels for each input cosmology.")
 63    else:
 64        Nreals = int(Nreals)
 65        Nthresh = int(Nthresh)
 66        b = np.stack(DensityFields)
 67        tmean = double([[0]*Nthresh]*b.shape[0])
 68        tstddev = double([[0]*Nthresh]*b.shape[0])
 69        n_cosmos, n_realizations, height, width = b.shape
 70        peaksarr = np.zeros((n_cosmos, Nthresh, Nreals), dtype=np.float64)
 71        thresh = np.linspace(-1, MaxThreshold, Nthresh)
 72        for cosmo in range(n_cosmos):
 73            for z in range(Nreals):
 74                bz = b[cosmo, z]
 75                local_max = (scipy.ndimage.maximum_filter(bz, size=3, mode='wrap')==bz)
 76                for t_idx, t in enumerate(thresh):
 77                    threshold_mask = (bz>=t)
 78                    peaks_mask = (local_max & threshold_mask)
 79                    peaksarr[cosmo, t_idx, z] = np.sum(peaks_mask)
 80            print('%d%% Done.'%(((double(cosmo)*double(Nreals))+double(z)+1.0)/(0.01*double(Nreals)*b.shape[0])))
 81        if(Type==0):
 82            for cosmo in range(0,b.shape[0],1):
 83                for t in range(0,Nthresh,1):
 84                    tmean[cosmo,t]=np.mean(peaksarr[cosmo,t])
 85                    tstddev[cosmo,t]=np.std(peaksarr[cosmo,t])
 86                if(Plot==1):
 87                    plt.plot(thresh,tmean[cosmo],label=CosmoLabels[cosmo])       
 88                    plt.fill_between(thresh, tmean[cosmo]-tstddev[cosmo], tmean[cosmo]+tstddev[cosmo], alpha=0.2)
 89                    plt.title('Peak Curves of Various Cosmologies')
 90        elif(Type==1):
 91            for cosmo in range(0,b.shape[0],1):
 92                for t in range(0,Nthresh,1):
 93                    tmean[cosmo,t]=np.mean(peaksarr[cosmo,t])/np.mean(peaksarr[0,t])
 94                    tstddev[cosmo,t]=np.std(peaksarr[cosmo,t])/np.mean(peaksarr[0,t])
 95                if(Plot==1):
 96                    plt.plot(thresh,tmean[cosmo],label=CosmoLabels[cosmo])       
 97                    plt.fill_between(thresh, tmean[cosmo]-tstddev[cosmo], tmean[cosmo]+tstddev[cosmo], alpha=0.2)
 98                    plt.title('Peak Curves of Various Cosmologies Normalized by ' + CosmoLabels[0])
 99        if(Plot==1):
100            plt.xlabel('Overdensity Threshold')
101            plt.ylabel('Number of Peaks')
102            if(LogScale==1):
103                plt.yscale('log')
104            plt.legend(bbox_to_anchor=(1.02,1),loc='upper left',borderaxespad=0)
105            plt.show()
106        return thresh, tmean, tstddev
107    
108    
109    
110def spherical_peaks(denslice,MaxThreshold,Nthresh):
111    # Peak Finder Helper Function
112    nside = int(np.sqrt(len(denslice)/12))
113    nops = 0
114    peaks = []
115    thresh = []
116    neigh = np.copy(hp.pixelfunc.get_all_neighbours(nside,range(len(denslice))))
117    for t in linspace(-1,MaxThreshold,Nthresh):
118        peax = np.copy(denslice[np.where((denslice[:]>=t)*(denslice[:]>=denslice[neigh[0,:]])*(denslice[:]>=denslice[neigh[1,:]])*(denslice[:]>=denslice[neigh[2,:]])*(denslice[:]>=denslice[neigh[3,:]])*(denslice[:]>=denslice[neigh[4,:]])*(denslice[:]>=denslice[neigh[5,:]])*(denslice[:]>=denslice[neigh[6,:]])*(denslice[:]>=denslice[neigh[7,:]]))])
119        nops = len(peax)
120        peaks.append(nops)
121        thresh.append(t)
122    return thresh,peaks
123
124def PeakCurves_Healpix(DensityFields=[], Nreals=10, MaxThreshold=16, Nthresh=101, Type=0, Plot=1, LogScale=1, CosmoLabels=['null']):
125    '''
126    Gives the peak curves for the given cosmologies' square (over)density fields.
127    
128    Parameters
129    ----------
130    DensityFields: float array of shape ''(nCosmo, NofRealizations, nPixels)''
131        3D Array of the healpix (over)density fields of the various cosmologies to be compared. The array should of shape (nCosmo, NofRealizations, nPixels) 
132        where 'nCosmo' is the number of cosmologies to be compared, 'NofRealizations' is the number of realizations of each input cosmology 
133        (NofRealizations >= Nreals), and 'nPixels' is the number of pixels in the healpix projected (over)density field such that [nPixels = 12*res*res], 
134        where 'res' is the resolution of the healpix map.
135        Example: np.array(DensityFields).shape = (3, 10, 512, 512) - 3 cosmologies containing 10 realisations each of (512x512) pixel 2D density fields.
136    Nreals: int
137        Number of realisations of the density fields to be used for the peak curves' statistics. Naturally, cannot be larger than the inherent number 
138        of realisations of the 2D density fields contained within the (.npy) data files input, i.e. in the above example, Nreals <=10. Any non-int values 
139        are type-cast to int.
140    MaxThreshold: float
141        Maximum overdensity threshold for which the peak values are to be plotted.
142    Nthresh: int
143        Number of overdensity threshold (X-axis) values to be computed in the closed interval [-1,MaxThreshold]. Any non-int values are type-cast to int.
144    Type: int: 0 or 1, optional
145        Type of peak curve to be plotted - 0 for raw peak curve plot; 1 for peak curves normalized by the first input cosmology's peak curve.
146    Plot: int: 0 or 1, optional
147        1 to plot the output and return results, 0 to skip plotting and only return results. Any values other than 1 will skip plotting.
148    LogScale: int: 0 or 1, optional
149        1 to plot the peak curves on log Y-axis, 0 to not use log Y-axis. Any values other than 1 will not output the log scale.
150    CosmoLabels: str array of shape ''(nCosmo)''
151        List of the names/labels to be assigned to the respective input cosmologies. Must have length equal to the number of cosmologies input ('nCosmo').
152        Example: ['CDM', 'SIDM', 'FDM']
153    
154    Returns
155    -------
156    thresh, tmean, tstddev: 3 numpy arrays containing the threshold values, mean peak values and std. dev. of the peak values
157        For (Type = 0): Peak curve plot of the various input comologies' density fields.
158        For (Type = 1): Peak curve plot of the various input comologies' density fields normalized by the first input cosmology's density field.
159        In both cases, the threshold values array (X-axis, 1D) 'thresh', the peaks array (Y-axis, 2D) 'tmean' containing the number of peaks corresponding
160        to the thresholds array for each input cosmology and their corresponding standard deviations (error, 2D) 'tstddev' are also returned.
161
162    Raises
163    ------
164    ValueError
165        If 'Type' is not 0 or 1.
166    ValueError
167        If 'Nreals' is lesser than 1.
168    ValueError
169        If 'MaxThreshold' is lesser than or equal to (-1).
170    ValueError
171        If the plot needs to be output and the number of labels ('CosmoLabels') is not equal to the number of input cosmologies.
172    '''
173    if((Type!=0) and (Type!=1)):
174        print("ERROR: Incorrect output type requested. Please use either 0 or 1 for variable 'Type'.")
175    elif(Nreals<=0):
176        print("ERROR: Invalid number of realizations entered. Please enter a positive integer value.")
177    elif(MaxThreshold<=-1):
178        print("ERROR: Please enter a maximum overdensity threshold value >(-1).")
179    elif((Plot==1) and (len(CosmoLabels)!=len(DensityFields))):
180        print("ERROR: Please enter the labels for each input cosmology.")
181    else:
182        Nreals = int(Nreals)
183        Nthresh = int(Nthresh)
184        b = np.stack(DensityFields)
185        peaksarr = double([[[0]*Nthresh]*Nreals]*b.shape[0])
186        tmean = double([[0]*Nthresh]*b.shape[0])
187        tstddev = double([[0]*Nthresh]*b.shape[0])
188        thresh = double([0]*Nthresh)
189        index = 0
190        for cosmo in range(0,b.shape[0],1):
191            for z in range(0,Nreals,1):
192                index += 1
193                thresh, peaksarr[cosmo,z] = spherical_peaks(b[cosmo,z],MaxThreshold,Nthresh)
194                if(index == b.shape[0]):
195                    print('%d%% Done.'%(((double(cosmo)*double(Nreals))+double(z)+1.0)/(0.01*double(Nreals)*b.shape[0])))
196                    index = 0
197        if(Type == 0):
198            for cosmo in range(0,b.shape[0],1):
199                tmean[cosmo] = np.mean(peaksarr[cosmo],axis=0)
200                tstddev[cosmo] = np.std(peaksarr[cosmo],axis=0)
201                if(Plot==1):
202                    plt.plot(thresh,tmean[cosmo],label=CosmoLabels[cosmo])       
203                    plt.fill_between(thresh, tmean[cosmo]-tstddev[cosmo], tmean[cosmo]+tstddev[cosmo], alpha=0.2)
204                    plt.title('Peak Curves of Various Cosmologies')
205        elif(Type == 1):
206            for cosmo in range(0,b.shape[0],1):
207                tmean[cosmo] = np.mean(peaksarr[cosmo],axis=0)/np.mean(peaksarr[0],axis=0)
208                tstddev[cosmo] = np.std(peaksarr[cosmo],axis=0)/np.mean(peaksarr[0],axis=0)
209                if(Plot==1):
210                    plt.plot(thresh,tmean[cosmo],label=CosmoLabels[cosmo])       
211                    plt.fill_between(thresh, tmean[cosmo]-tstddev[cosmo], tmean[cosmo]+tstddev[cosmo], alpha=0.2)
212                    plt.title('Peak Curves of Various Cosmologies Normalized by ' + CosmoLabels[0])
213        if(Plot==1):
214            plt.xlabel('Overdensity Threshold')
215            plt.ylabel('Number of Peaks')
216            if(LogScale == 1):
217                plt.yscale('log')
218            plt.legend(bbox_to_anchor=(1.02,1),loc='upper left',borderaxespad=0)
219            plt.show()
220        return thresh, tmean, tstddev