__author__ = 'JB'
import warnings
from scipy.optimize import leastsq
from scipy.interpolate import interp1d
from astropy.modeling import models, fitting
from pyklip.kpp.utils.mathfunc import *
from pyklip.kpp.utils.GPIimage import *
[docs]
def get_image_stat_map(image,
image_without_planet=None,
IOWA = None,
N = None,
centroid = None,
r_step = 2,
Dr = 2,
type = "SNR",
image_wide = None,
azimuthal_sectors_for_masking=1):
"""
Calculate the SNR, the standard deviation or the probability (tail distribution) of a given image using concentric
annuli.
Args:
image: The image for which one wants the statistic.
image_without_planet: Same as image but where any real signal has been masked out. The code will use
"image_without_planet" to calculate the standard deviation or the PDF.
This can be a negatively derotated image in the context of ADI.
IOWA: (IWA,OWA) inner working angle, outer working angle. It defines boundary to the zones in which the
statistic is calculated.
If None, kpp.utils.GPIimage.get_IOWA() is used.
N: Defines the width of the ring by fixing the number of pixels of the annulus.
The width of the annuli will therefore vary with sepration.
centroid: (x_cen,y_cen) Define the center of the image.
Default is x_cen = (nx-1)//2 ; y_cen = (ny-1)//2
r_step: (default=2pix) Distance between two consecutive annuli mean separation (in pixel).
Dr: (default=2pix) Width of the annulus (in pixel).
type: Indicate the type of statistic to be calculated.
If "SNR" (default) simple stddev calculation and returns SNR.
If "stddev" returns the pure standard deviation map.
If "mean" returns a map from the radial mean.
If "sum" returns a map from the radial sum.
If "proba" triggers proba calculation with pdf fitting.
image_wide: Don't divide the image in annuli or sectors when computing the statistic.
Use the entire image directly.
azimuthal_sectors_for_masking: if > 1, breaks up the image into this many azimuthal
sectors, with pixels in the sector and half a sector outside all
masked to compute the statistic in that sector. This is to avoid planet
signals from downweighter their own SNR. Useful when no image without
planet is available initially. Doesn't work for proba.
Return:
The statistic map for image.
"""
if image_without_planet is None:
image_without_planet = image
ny,nx = image.shape
if centroid is None :
centroid = ((nx-1)//2 ,(ny-1)//2)
if image_wide is None:
image_wide = False
if IOWA is None:
IWA,OWA = get_IOWA(image_without_planet, centroid = centroid)
else:
IWA,OWA = IOWA
if type == "proba":
pdf_list, cdf_list, sampling_list, annulus_radii_list = get_image_PDF(image_without_planet,(IWA,OWA),N,centroid,r_step=r_step,Dr=Dr,image_wide = image_wide)
pdf_radii = np.array(annulus_radii_list)[:,0]
stat_map = np.zeros(image.shape) + np.nan
# Build the x and y coordinates grids
x_grid, y_grid = np.meshgrid(np.arange(nx)-centroid[0], np.arange(ny)-centroid[1])
# Calculate the radial distance of each pixel
r_grid = abs(x_grid +y_grid*1j)
image_finite = np.where(np.isfinite(image))
#Build the cdf_models from interpolation
cdf_interp_list = []
for sampling,cdf_sampled in zip(sampling_list,cdf_list):
cdf_interp_list.append(interp1d(sampling,cdf_sampled,kind = "linear",bounds_error = False, fill_value=1.0))
#f = interp1d(sampling,cdf_sampled,kind = "linear",bounds_error = False, fill_value=1.0)
#plt.plot(np.arange(-10,10,0.1),f(np.arange(-10,10,0.1)))
#plt.show()
for k,l in zip(image_finite[0],image_finite[1]):
#stdout.flush()
#stdout.write("\r%d" % k)
r = r_grid[k,l]
if r < OWA:
r_closest_id, r_closest = min(enumerate(pdf_radii), key=lambda x: abs(x[1]-r))
if (r-r_closest) < 0:
r_closest_id2 = r_closest_id - 1
else:
r_closest_id2 = r_closest_id + 1
if (r_closest_id2 < 0) or (r_closest_id2 > (pdf_radii.size-1)):
stat_map[k,l] = 1-cdf_interp_list[r_closest_id](image[k,l])
#plt.plot(np.arange(-10,10,0.1),cdf(np.arange(-10,10,0.1)))
#plt.show()
else:
r_closest2 = pdf_radii[r_closest_id2]
stat_map[k,l] = 1-(cdf_interp_list[r_closest_id](image[k,l])*abs(r-r_closest2)+cdf_interp_list[r_closest_id2](image[k,l])*abs(r-r_closest))/abs(r_closest-r_closest2)
else:
stat_map[k,l] = 1-cdf_interp_list[pdf_radii.size-1](image[k,l])
if 0:
import matplotlib.pyplot as plt
plt.figure(1)
plt.subplot(1,3,1)
plt.imshow(np.log10(stat_map),interpolation="nearest")
plt.colorbar()
plt.subplot(1,3,2)
plt.imshow(image,interpolation="nearest")
plt.subplot(1,3,3)
plt.imshow(image_without_planet,interpolation="nearest")
plt.show()
return -np.log10(stat_map)
else:
if type =="SNR":
tmp_type = "stddev"
else:
tmp_type = type
azimuthal_bounds = np.linspace(0, 360, azimuthal_sectors_for_masking+1)
stat_map = np.zeros(image.shape) + np.nan
ny,nx = image.shape
# Build the x and y coordinates grids
x_grid, y_grid = np.meshgrid(np.arange(nx)-centroid[0], np.arange(ny)-centroid[1])
# Calculate the radial distance of each pixel
r_grid = abs(x_grid +y_grid*1j)
theta_grid = np.degrees(np.arctan2(y_grid, x_grid)) % 360
for azimuthal_index in range(azimuthal_sectors_for_masking):
theta_start = azimuthal_bounds[azimuthal_index]
theta_end = azimuthal_bounds[azimuthal_index+1]
dtheta = theta_end - theta_start
in_zone = np.where((theta_grid >= theta_start) & (theta_grid < theta_end))
if azimuthal_sectors_for_masking > 1:
image_without_planet_azimuthal_mask = np.copy(image_without_planet)
azimuthal_mask = np.where((theta_grid >= theta_start - dtheta/2) & (theta_grid < theta_end + dtheta/2))
image_without_planet_azimuthal_mask[azimuthal_mask] = np.nan
if theta_start == 0:
image_without_planet_azimuthal_mask[np.where(theta_grid >= 360 - dtheta/2)] = np.nan
if theta_end == 360:
image_without_planet_azimuthal_mask[np.where(theta_grid < dtheta/2)] = np.nan
else:
image_without_planet_azimuthal_mask = image_without_planet
stat_list, annulus_radii_list = get_image_stat(image_without_planet_azimuthal_mask,tmp_type,(IWA,OWA),N,centroid,r_step=r_step,Dr=Dr,
image_wide=image_wide)
radii = np.array(annulus_radii_list)[:,0]
#print(radii,stddev_list)
if not image_wide:
stat_func = interp1d(radii,stat_list,kind = "linear",bounds_error = False, fill_value=np.nan)
else:
stat_func = lambda x: stat_list[0]
#plt.figure()
#plt.plot(np.linspace(0,140,200),stddev_func(np.linspace(0,140,200)))
#plt.show()
stat_vals = stat_func(r_grid[in_zone].ravel())
stat_map[in_zone] = stat_vals
# stat_map = stat_func(r_grid.ravel())
# stat_map.shape = r_grid.shape
nanpix = np.where(np.isnan(image))
stat_map[nanpix] = np.nan
if type == "SNR":
stat_map = image/stat_map
return stat_map
image_finite = np.where(np.isfinite(image))
for k,l in zip(image_finite[0],image_finite[1]):
#stdout.flush()
#stdout.write("\r%d" % k)
r = r_grid[k,l]
if type == "SNR":
stat_map[k,l] = image[k,l]/stddev_func(r)
elif type == "stddev":
stat_map[k,l] = stddev_func(r)
return stat_map
[docs]
def get_image_PDF(image,IOWA=None,N = 2000,centroid = None, r_step = None,Dr=None,image_wide = None):
"""
Calculate the PDF of a given image using annuli.
Args:
image: The image for which one wants the statistic.
IOWA: (IWA,OWA) inner working angle, outer working angle. It defines boundary to the zones in which the
statistic is calculated.
If None, kpp.utils.GPIimage.get_IOWA() is used.
N: Defines the width of the ring by the number of pixels it has to include.
The width of the annuli will therefore vary with sepration. Default is N=3000.
centroid: (x_cen,y_cen) Define the center of the image.
Default is x_cen = (nx-1)//2 ; y_cen = (ny-1)//2
r_step: Distance between two consecutive annuli mean separation (in pixel).
Dr: Width of the annulus (in pixel).
image_wide: Don't divide the image in annuli or sectors when computing the statistic.
Use the entire image directly. Not available if "pixel based: is defined,
Return:
pdf_list: List of PDF values for each annulus. The sampling of each PDF can be found in sampling_list.
cdf_list: CDF values for each annulus. The sampling of each CDF can be found in sampling_list.
sampling_list: Sampling for the PDF and the CDF.
annulus_radii_list: List of ((r_min+r_max)/2.,r_min,r_max) with r_min,r_max the boundaries of an annulus.
"""
if image_wide is None:
image_wide = False
if IOWA is None:
IWA,OWA = get_IOWA(image, centroid = centroid)
else:
IWA,OWA = IOWA
ny,nx = image.shape
if 0:
import matplotlib.pyplot as plt
fig = 1
plt.figure(fig,figsize=(16,8))
plt.subplot(121)
plt.imshow(image,interpolation="nearest")
plt.colorbar()
data = image[np.where(np.isfinite(image))]
im_std = np.std(data)
bins = np.arange(np.min(data),np.max(data),im_std/10.)
im_histo = np.histogram(data, bins=bins)[0]
N_bins = bins.size-1
center_bins = 0.5*(bins[0:N_bins]+bins[1:N_bins+1])
plt.subplot(122)
plt.plot(center_bins,np.array(im_histo,dtype="double"),'bx-', markersize=5,linewidth=3)
plt.grid(True)
ax = plt.gca()
ax.set_yscale('log')
plt.show()
image_mask = np.ones((ny,nx))
image_mask[np.where(np.isnan(image))] = 0
if centroid is None :
x_cen = (nx-1)//2 ; y_cen = (ny-1)//2
else:
x_cen, y_cen = centroid
# Build the x and y coordinates grids
x, y = np.meshgrid(np.arange(nx)-x_cen, np.arange(ny)-y_cen)
# Calculate the radial distance of each pixel
r_grid = abs(x +y*1j)
#th_grid = np.arctan2(x,y)
if not image_wide:
# Define the radii intervals for each annulus
if Dr is None:
r0 = IWA
annuli_radii = []
if r_step is None:
while np.sqrt(N/np.pi+r0**2) < OWA:
annuli_radii.append((r0,np.sqrt(N/np.pi+r0**2)))
r0 = np.sqrt(N/np.pi+r0**2)
else:
while np.sqrt(N/np.pi+r0**2) < OWA:
annuli_radii.append((r0,np.sqrt(N/np.pi+r0**2)))
r0 += r_step
annuli_radii.append((r0,np.max([ny,nx])))
else:
annuli_radii = []
for r in np.arange(IWA,OWA+Dr,r_step):
annuli_radii.append((r-Dr/2.,r+Dr/2.))
else:
annuli_radii = [(IWA,OWA)]
N_annuli = len(annuli_radii)
pdf_list = []
cdf_list = []
sampling_list = []
annulus_radii_list = []
if 0:
rings = np.zeros((ny,nx))+np.nan
for it, rminmax in enumerate(annuli_radii):
r_min,r_max = rminmax
#print(rminmax)
where_ring = np.where((r_min< r_grid) * (r_grid < r_max) * image_mask)
#print(np.size(where_ring[0]))
if 0:
import matplotlib.pyplot as plt
image_tmp = copy(image)
image_tmp[where_ring] = np.nan
plt.figure(2)
plt.imshow(image_tmp,interpolation="nearest")
plt.show()
if 0:
rings[where_ring] = it
data = image[where_ring]
cdf_model, pdf_model, sampling, im_histo, center_bins = get_cdf_model(data)
pdf_list.append(pdf_model)
cdf_list.append(cdf_model)
sampling_list.append(sampling)
annulus_radii_list.append(((r_min+r_max)/2.,r_min,r_max))
if 0:
import matplotlib.pyplot as plt
fig = 1
plt.figure(fig,figsize=(8,8))
plt.subplot(np.ceil(np.sqrt(N_annuli)),np.ceil(np.sqrt(N_annuli)),it)
plt.plot(sampling,pdf_model,'b-',linewidth=3)
plt.plot(sampling,1.-cdf_model,'r-',linewidth=3)
plt.xlabel('criterion value', fontsize=20)
plt.ylabel('Probability of the value', fontsize=20)
plt.grid(True)
ax = plt.gca()
ax.tick_params(axis='x', labelsize=20)
ax.tick_params(axis='y', labelsize=20)
ax.legend(['flat cube histogram','flat cube histogram (Gaussian fit)','planets'], loc = 'upper right', fontsize=12)
ax.set_yscale('log')
plt.ylim((10**-7,10))
if 0:
import matplotlib.pyplot as plt
plt.figure(2,figsize=(8,8))
plt.imshow(rings,interpolation="nearest")
plt.show()
return pdf_list, cdf_list, sampling_list, annulus_radii_list
[docs]
def get_image_stddev(image,
IOWA = None,
N = None,
centroid = None,
r_step = 2,
Dr=2,
image_wide = None):
"""
Calculate the standard deviation of a given image using annuli.
Args:
image: The image for which one wants the statistic.
IOWA: (IWA,OWA) inner working angle, outer working angle. It defines boundary to the zones in which the
statistic is calculated.
If None, kpp.utils.GPIimage.get_IOWA() is used.
N: Defines the width of the ring by the number of pixels it has to include.
The width of the annuli will therefore vary with sepration. Default is N=3000.
centroid: (x_cen,y_cen) Define the center of the image.
Default is x_cen = (nx-1)//2 ; y_cen = (ny-1)//2
r_step: (default=2pix) Distance between two consecutive annuli mean separation (in pixel).
Dr: (default=2pix) Width of the annulus (in pixel).
image_wide: Don't divide the image in annuli or sectors when computing the statistic.
Use the entire image directly. Not available if "pixel based: is defined,
Return:
stddev_list: standard deviation values at the center of each
annulus_radii_list: List of ((r_min+r_max)/2.,r_min,r_max) with r_min,r_max the boundaries of an annulus.
"""
return get_image_stat(image,"stddev",
IOWA = IOWA,
N = N,
centroid = centroid,
r_step = r_step,
Dr=Dr,
image_wide = image_wide)
[docs]
def get_image_stat(image,type,
IOWA = None,
N = None,
centroid = None,
r_step = 2,
Dr=2,
image_wide = None):
"""
Calculate the standard deviation or the mean of a given image using annuli.
Args:
image: The image or cubes for which one wants the statistic.
type: "stddev" or "mean" or "sum"
IOWA: (IWA,OWA) inner working angle, outer working angle. It defines boundary to the zones in which the
statistic is calculated.
If None, kpp.utils.GPIimage.get_IOWA() is used.
N: Defines the width of the ring by the number of pixels it has to include.
The width of the annuli will therefore vary with sepration. Default is N=3000.
centroid: (x_cen,y_cen) Define the center of the image.
Default is x_cen = (nx-1)//2 ; y_cen = (ny-1)//2
r_step: Distance between two consecutive annuli mean separation. Not available if "pixel based" is defined,
Dr: If not None defines the width of the ring as Dr. N is then ignored if Dth is defined.
image_wide: Don't divide the image in annuli or sectors when computing the statistic.
Use the entire image directly. Not available if "pixel based: is defined,
Return:
stddev_list: standard deviation values at the center of each
annulus_radii_list: List of ((r_min+r_max)/2.,r_min,r_max) with r_min,r_max the boundaries of an annulus.
"""
if image_wide is None:
image_wide = False
if IOWA is None:
IWA,OWA,inner_mask,outer_mask = get_occ(image, centroid = centroid)
else:
IWA,OWA = IOWA
ny,nx = image.shape
image_mask = np.ones((ny,nx))
image_mask[np.where(np.isnan(image))] = 0
if centroid is None :
x_cen = (nx-1)//2 ; y_cen = (ny-1)//2
else:
x_cen, y_cen = centroid
# Build the x and y coordinates grids
x, y = np.meshgrid(np.arange(nx)-x_cen, np.arange(ny)-y_cen)
# Calculate the radial distance of each pixel
r_grid = abs(x +y*1j)
#th_grid = np.arctan2(x,y)
if not image_wide:
# Define the radii intervals for each annulus
if Dr is None:
r0 = IWA
annuli_radii = []
if r_step is None:
while np.sqrt(N/np.pi+r0**2) < OWA:
annuli_radii.append((r0,np.sqrt(N/np.pi+r0**2)))
r0 = np.sqrt(N/np.pi+r0**2)
else:
while np.sqrt(N/np.pi+r0**2) < OWA:
annuli_radii.append((r0,np.sqrt(N/np.pi+r0**2)))
r0 += r_step
annuli_radii.append((r0,np.max([ny,nx])))
else:
annuli_radii = []
for r in np.arange(IWA,OWA+Dr,r_step):
annuli_radii.append((r-Dr/2.,r+Dr/2.))
else:
annuli_radii = [(IWA,OWA)]
#N_annuli = len(annuli_radii)
stat_list = []
annulus_radii_list = []
for it, rminmax in enumerate(annuli_radii):
r_min,r_max = rminmax
where_ring = np.where((r_min< r_grid) * (r_grid < r_max) * image_mask)
data = image[where_ring]
if type == "stddev":
stat = np.nanstd(data)
elif type == "mean":
stat = np.nanmean(data)
elif type == "sum":
stat = np.nansum(data)
stat_list.append(stat)
annulus_radii_list.append(((r_min+r_max)/2.,r_min,r_max))
return stat_list, annulus_radii_list
[docs]
def get_cdf_model(data,interupt_plot = False,pure_gauss=False):
"""
Calculate a model CDF for some data.
/!\ This function is for some reason still a work in progress. JB could never decide what the best option was.
But it should work even if the code is a mess.
Args:
data: arrays of samples from a random variable
interupt_plot: Plot the histogram and model fit. It
pure_gauss: Assume gaussian statistic. Do not fit exponential tails.
Return: (cdf_model,new_sampling,im_histo, center_bins) with:
cdf_model: The cdf model = np.cumsum(pdf_model)
pdf_model: The pdf model
sampling: sampling of pdf/cdf_model
im_histo: histogram from original data
center_bins: bin centers for im_histo
"""
pdf_model,sampling,im_histo,center_bins = get_pdf_model(data,interupt_plot=interupt_plot,pure_gauss=pure_gauss)
return np.cumsum(pdf_model),pdf_model,sampling,im_histo,center_bins
[docs]
def get_pdf_model(data,interupt_plot = False,pure_gauss = False):
"""
Calculate a model PDF for some data.
/!\ This function is for some reason still a work in progress. JB could never decide what the best option was.
But it should work even if the code is a mess.
Args:
data: arrays of samples from a random variable
interupt_plot: Plot the histogram and model fit. It
pure_gauss: Assume gaussian statistic. Do not fit exponential tails.
Return: (pdf_model,new_sampling,im_histo, center_bins) with:
pdf_model: The pdf model
new_sampling: sampling of pdf_model
im_histo: histogram from original data
center_bins: bin centers for im_histo
"""
im_std = np.std(data)
#print(im_std)
bins = np.arange(np.min(data),np.max(data),im_std/5.)
im_histo = np.histogram(data, bins=bins)[0]
N_bins = bins.size-1
center_bins = 0.5*(bins[0:N_bins]+bins[1:N_bins+1])
g_init = models.Gaussian1D(amplitude=np.max(im_histo), mean=0.0, stddev=im_std)
fit_g = fitting.LevMarLSQFitter()
warnings.simplefilter('ignore')
g = fit_g(g_init, center_bins, im_histo)#, weights=1/im_histo)
g.stddev = abs(g.stddev)
right_side_noZeros = np.where((center_bins > (g.mean+2*g.stddev))*(im_histo != 0))
N_right_bins_noZeros = len(right_side_noZeros[0])
left_side_noZeros = np.where((center_bins < (g.mean-2*g.stddev))*(im_histo != 0))
N_left_bins_noZeros = len(left_side_noZeros[0])
right_side = np.where((center_bins > (g.mean+2*g.stddev)))
left_side = np.where((center_bins < (g.mean-2*g.stddev)))
if not pure_gauss:
if N_right_bins_noZeros < 5:
where_pos_zero = np.where((im_histo == 0) * (center_bins > g.mean))
if len(where_pos_zero[0]) != 0:
right_side_noZeros = (range(where_pos_zero[0][0]-5,where_pos_zero[0][0]),)
right_side = (range(where_pos_zero[0][0]-5,center_bins.size),)
else:
right_side_noZeros = (range(center_bins.size-5,center_bins.size),)
right_side = right_side_noZeros
N_right_bins_noZeros = 5
if N_left_bins_noZeros < 5:
where_neg_zero = np.where((im_histo == 0) * (center_bins < g.mean))
if len(where_neg_zero[0]) != 0:
left_side_noZeros = (range(where_neg_zero[0][len(where_neg_zero[0])-1]+1,where_neg_zero[0][len(where_neg_zero[0])-1]+6),)
left_side = (range(0,where_neg_zero[0][len(where_neg_zero[0])-1]+6),)
else:
left_side_noZeros = (range(0,5),)
left_side = left_side_noZeros
N_left_bins_noZeros = 5
#print(left_side,right_side)
#print(im_histo[left_side],im_histo[right_side])
#print(right_side_noZeros,left_side_noZeros)
#print(im_histo[right_side_noZeros],im_histo[left_side_noZeros])
#print(N_right_bins_noZeros,N_left_bins_noZeros)
if N_right_bins_noZeros >= 2:
alpha0 = (np.log(im_histo[right_side_noZeros[0][N_right_bins_noZeros-1]])-np.log(im_histo[right_side_noZeros[0][0]]))/(center_bins[right_side_noZeros[0][0]]-center_bins[right_side_noZeros[0][N_right_bins_noZeros-1]])
m_alpha0 = -np.log(im_histo[right_side_noZeros[0][0]])-alpha0*center_bins[right_side_noZeros[0][0]]
param0_rightExp = (m_alpha0,alpha0)
LSQ_func = lambda para: LSQ_model_exp((bins[0:bins.size-1])[right_side], im_histo[right_side],para[0],para[1])
param_fit_rightExp = leastsq(LSQ_func,param0_rightExp)
else:
param_fit_rightExp = None
#print(param0_rightExp,param_fit_rightExp)
if N_left_bins_noZeros >= 2:
alpha0 = (np.log(im_histo[left_side_noZeros[0][N_left_bins_noZeros-1]])-np.log(im_histo[left_side_noZeros[0][0]]))/(center_bins[left_side_noZeros[0][0]]-center_bins[left_side_noZeros[0][N_left_bins_noZeros-1]])
m_alpha0 = -np.log(im_histo[left_side_noZeros[0][0]])-alpha0*center_bins[left_side_noZeros[0][0]]
param0_leftExp = (m_alpha0,alpha0)
LSQ_func = lambda para: LSQ_model_exp((bins[0:bins.size-1])[left_side], im_histo[left_side],para[0],para[1])
param_fit_leftExp = leastsq(LSQ_func,param0_leftExp)
else:
param_fit_leftExp = None
#print(param0_leftExp,param_fit_leftExp)
new_sampling = np.arange(2*np.min(data),4*np.max(data),im_std/100.)
if pure_gauss:
pdf_model = g(new_sampling)
pdf_model_exp = new_sampling*0
else:
pdf_model_gaussian = interp1d(center_bins,np.array(im_histo,dtype="double"),kind = "cubic",bounds_error = False, fill_value=0.0)(new_sampling)
if not pure_gauss:
right_side2 = np.where((new_sampling >= g.mean))
left_side2 = np.where((new_sampling < g.mean))
#print(g.mean+0.0,g.stddev+0.0)
pdf_model_exp = np.zeros(new_sampling.size)
weights = np.zeros(new_sampling.size)
if param_fit_rightExp is not None:
pdf_model_exp[right_side2] = model_exp(new_sampling[right_side2],*param_fit_rightExp[0])
weights[right_side2] = np.tanh((new_sampling[right_side2]-(g.mean+2*g.stddev))/(0.1*g.stddev))
else:
weights[right_side2] = -1.
if param_fit_leftExp is not None:
pdf_model_exp[left_side2] = model_exp(new_sampling[left_side2],*param_fit_leftExp[0])
weights[left_side2] = np.tanh(-(new_sampling[left_side2]-(g.mean-2*g.stddev))/(0.1*g.stddev))
else:
weights[left_side2] = -1.
weights = 0.5*(weights+1.0)
#weights[np.where(weights > 1-10^-3)] = 1
pdf_model = weights*pdf_model_exp + (1-weights)*pdf_model_gaussian
#pdf_model[np.where(weights > 1-10^-5)] = pdf_model_exp[np.where(pdf_model > 1-10^-5)]
if 0:
import matplotlib.pyplot as plt
fig = 2
plt.figure(fig,figsize=(8,8))
plt.plot(new_sampling, weights, "r")
#plt.plot(new_sampling, (1-weights), "--r")
#plt.plot(new_sampling, pdf_model_exp, "g")
#plt.plot(new_sampling, pdf_model_gaussian, "b")
#plt.plot(new_sampling, pdf_model, "c") #/np.sum(pdf_model)
#plt.plot(new_sampling, 1-np.cumsum(pdf_model/np.sum(pdf_model)), "--.")
ax = plt.gca()
#ax.set_yscale('log')
plt.grid(True)
#plt.ylim((10**-15,100000))
#plt.xlim((1*np.min(data),2*np.max(data)))
plt.show()
if interupt_plot:
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams.update({'font.size': 20})
fig = 2
plt.close(2)
plt.figure(fig,figsize=(16,8))
plt.subplot(121)
plt.plot(new_sampling,pdf_model,'r-',linewidth=5)
plt.plot(center_bins,g(center_bins),'c--',linewidth=3)
plt.plot(new_sampling,pdf_model_exp,'g--',linewidth=3)
plt.plot(center_bins,np.array(im_histo,dtype="double"),'b.', markersize=10,linewidth=3)
#plt.plot(new_sampling,np.cumsum(pdf_model),'g.')
plt.xlabel('Metric value')
plt.ylabel('Number per bin')
plt.xlim((2*np.min(data),2*np.max(data)))
plt.grid(True)
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
ax = plt.gca()
ax.tick_params(axis='x')
ax.tick_params(axis='y')
ax.legend(['PDF Model Fit','Central Gaussian Fit','Tails Exponential Fit','Histogram'], loc = 'lower left', fontsize=15)
ax.set_yscale('log')
plt.ylim((10**-1,10000))
pdf_model /= np.sum(pdf_model)
if interupt_plot:
from mpl_toolkits.axes_grid1 import host_subplot
import mpl_toolkits.axisartist as AA
host = host_subplot(122, axes_class=AA.Axes)
par1 = host.twinx()
p1, = host.plot(new_sampling,pdf_model/(new_sampling[1]-new_sampling[0]),'r-',linewidth=5)
host.tick_params(axis='x', labelsize=20)
host.tick_params(axis='y', labelsize=20)
host.set_ylim((10**-3,10**2))
host.set_yscale('log')
p2, = par1.plot(new_sampling,1-np.cumsum(pdf_model),'g-',linewidth=5)
par1.set_ylabel("False positive rate")
par1.set_yscale('log')
par1.set_ylim((10**-4,10.))
host.axis["left"].label.set_color(p1.get_color())
par1.axis["right"].label.set_color(p2.get_color())
plt.xlabel('Metric value')
plt.ylabel('Probability density')
plt.xlim((2*np.min(data),2*np.max(data)))
plt.grid(True)
plt.legend(['PDF model','Tail distribution'], loc = 'lower left', fontsize=15)
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.show()
return pdf_model,new_sampling,np.array(im_histo,dtype="double"), center_bins
[docs]
def get_cube_stddev(cube,IOWA,N = 2000,centroid = None, r_step = None,Dr=None):
# Not tested
nl,ny,nx = cube.shape
stddev_table = []
annulus_radii_table = []
for k in range(nl):
stddev_list, annulus_radii_list = get_image_stddev(cube[k,:,:],IOWA,N = N,centroid = centroid, r_step = r_step,Dr=Dr)
stddev_table.append(stddev_list)
annulus_radii_table.append(annulus_radii_list)
return stddev_table,annulus_radii_table