Commit 02bb0e88 authored by LUSTIG Peter's avatar LUSTIG Peter

Purity and completness (completness not testet yet) for masked region

parent a239436a
......@@ -27,16 +27,21 @@ from matplotlib.ticker import FormatStrFormatter
from collections import OrderedDict
from utils import completness_purity_wcs, completness_worker, purity_worker
from utils import find_nearest
from nikamap.utils import pos_in_mask
class PCEvaluation:
def __init__(self, sources, fake_sources=None, shape=None, wcs=None,
def __init__(self, sources, fake_sources=None, realsources=None,
shape=None, wcs=None,
flux=None, mapbins=19, threshold_bins=5,
threshold_range=(3, 5)):
self.flux = flux
self.sources = sources
self.fake_sources = fake_sources
self.real_sources = realsources
self.real_sources_map = None
self.threshold_range = threshold_range
if self.flux is not None:
idxsort = np.argsort(flux.to_value(u.mJy))
......@@ -69,6 +74,9 @@ class PCEvaluation:
bins=self.bins,
threshold_range=threshold_range,
threshold_bins=threshold_bins)
self.shape3D = (mapbins, mapbins, threshold_bins)
print("shape", self.shape_3D)
print("wcs", self.wcs_3D)
# Testing the lower edges
wcs_threshold = self.wcs_3D.sub([3])
......@@ -80,6 +88,10 @@ class PCEvaluation:
self.purity,
self.hitmap,
self.detectmap) = self.GetCP()
(_, _, _,
self.real_sources_map) = self.completness_purity(
self.real_sources, None,
self.wcs_3D, self.shape_3D)
def GetCompletnessBin(self, xbin, ybin):
return self.completness[:, ybin, xbin, :]
......@@ -291,14 +303,214 @@ class PCEvaluation:
data = data[:, :, threshidx]
thresh = self.thresholds[threshidx]
data /= normalize
print(data.shape)
plt.figure()
plt.title('Detected Sources (SNR {})'.format(thresh), fontsize=30)
plt.imshow(data, origin='lower', aspect='auto', **kwargs)
plt.imshow(data/normalize, origin='lower', **kwargs)
cbar = plt.colorbar()
cbar.ax.tick_params(labelsize=tickfs)
def PlotDarkcountsThresholds(self, xbin, ybin, fluxbin=None, normalize=1):
data = self.detectmap[:, ybin, xbin, :]
if fluxbin is None:
data = np.sum(data, axis=0)
ax = Plot1D(self.thresholds, data/normalize, title='Absolute Counts',
xlabel='Detection Threshold',
ylabel='Counts', label='Jackknife')
if self.real_sources_map is not None:
rdata = self.real_sources_map[ybin, xbin, :]
ax.plot(self.thresholds, rdata, label='Map')
ax.legend(loc='upper right', fontsize=20)
ax.set_yscale("log", nonposy='clip')
class PCAreaEvaluation:
def __init__(self, mask=None, sources=None, fake_sources=None, wcs=None,
threshold_bins=5, threshold_range=(3, 5), realsources=None,
nsim=None):
assert sources is not None, 'Source catalog must be provided'
self.mask = mask
self.threshold_range = threshold_range
self.threshold_bins = threshold_bins
self.thresholds = np.linspace(threshold_range[0], threshold_range[1],
threshold_bins)
self.nsim = nsim
self.real_sources = realsources
sources = AddXYCoordinatesTo(sources, wcs)
realsources = AddXYCoordinatesTo(realsources, wcs)
fake_sources = AddXYCoordinatesTo(fake_sources, wcs)
if mask is not None:
assert wcs is not None, 'Need wcs to find sources in mask.'
self.sources = GetInsideSources(sources, self.mask)
if realsources is not None:
realsources = GetInsideSources(realsources, self.mask)
self.RealSources = realsources
_, self.absolute_detections = Purity1D(self.sources,
self.threshold_range,
self.threshold_bins)
self.normed_detections = self.GetNormedDetections(
self.absolute_detections,
self.nsim)
_, self.RealmapDetections = Purity1D(self.RealSources,
self.threshold_range,
self.threshold_bins)
def GetNormedDetections(self, absolute_detections=None, norm=None):
if absolute_detections is None:
absolute_detections = self.absolute_detections
if absolute_detections is not None and norm is not None:
return absolute_detections / norm
else:
return None
def PlotMask(self, mask=None, sources=False, sourcecat=None,
realsources=False, realsourcecat=None, **kwargs):
if mask is None:
mask = self.mask
if sources and sourcecat is None:
sourcecat = self.sources
if realsources and realsourcecat is None:
sourcecat = self.RealSources
ax, img = Plot2D(np.array(mask, dtype=int), **kwargs)
if sourcecat is not None:
plt.scatter(sourcecat['xmean'], sourcecat['ymean'], marker='.')
if realsourcecat is not None:
plt.scatter(realsourcecat['xmean'], realsourcecat['ymean'],
marker='.', color='r')
return ax, img
def PlotDarkcountsThresholds(self, data=None, fluxbin=None, normalize=1):
if data is None:
if self.normed_detections is not None:
data = self.normed_detections
else:
data = self.absolute_detections
print(data.shape)
print(self.thresholds.shape)
if fluxbin is None and len(data.shape) == 2:
data = np.sum(data, axis=0)
elif fluxbin is not None and len(data.shape) == 2:
data = data[fluxbin]
print(data.shape)
ax = Plot1D(self.thresholds, data/normalize, title='Absolute Counts',
xlabel='Detection Threshold',
ylabel='Counts', label='Jackknife')
if self.RealmapDetections is not None:
ax.plot(self.thresholds, self.RealmapDetections, label='Map')
ax.legend(loc='upper right', fontsize=20)
'''
plt.figure()
plt.plot(self.thresholds, 1./(self.RealmapDetections / (data/normalize)))
plt.show(block=True)
'''
ax.set_yscale("log", nonposy='clip')
def Completness1D(sources, fake_sources, threshold_range, nbins):
if fake_sources is not None:
min_threshold, max_threshold = threshold_range
fake_snr = np.ma.array(sources[fake_sources['find_peak']]['SNR'],
mask=fake_sources['find_peak'].mask)
fake_snr[fake_snr > max_threshold] = max_threshold
bin_edges = np.linspace(min_threshold, max_threshold, nbins)
# add right edge to last bin
bin_edges = np.append(bin_edges, bin_edges[-1]+1)
completness, _ = np.histogram(fake_snr, bins=bin_edges,
weights=~fake_sources['find_peak'].mask)
completness = np.cumsum(completness[::-1])[::-1]
norm_comp = len(fake_sources)
else:
completness, norm_comp = None, None
return completness, norm_comp
def Purity1D(sources, threshold_range, nbins):
if sources is not None:
min_threshold, max_threshold = threshold_range
sources_snr = sources['SNR']
sources_snr[sources_snr > max_threshold] = max_threshold
bin_edges = np.linspace(min_threshold, max_threshold, nbins)
# add right edge to last bin
bin_edges = np.append(bin_edges, bin_edges[-1]+1)
if 'fake_sources' in sources.keys():
purity, _ = np.histogram(sources_snr, bins=bin_edges,
weights=~sources['fake_sources'].mask)
purity = np.cumsum(purity[::-1])[::-1]
else:
purity = None
norm_pur, _ = np.histogram(sources_snr, bins=bin_edges)
norm_pur = np.cumsum(norm_pur[::-1])[::-1]
else:
purity, norm_pur = None, None
return purity, norm_pur
def AddXYCoordinatesTo(sources, wcs):
if sources is not None:
pos = np.array(wcs.wcs_world2pix(sources['ra'],
sources['dec'], 0)).T
sources['xmean'], sources['ymean'] = pos[:, 0], pos[:, 1]
return sources
def GetInsideSources(sources, mask):
pos = np.array([sources['xmean'], sources['ymean']]).T
pos, sources = pos_in_mask(pos, mask=mask, extra=sources)
return sources
def Plot1D(x, y, title='', xlabel='', ylabel='', **kwargs):
ax = plt.subplot(111)
ax.plot(x, y, **kwargs)
ax.set_title(title, fontsize=30, y=1.02)
ax.set_xlabel(xlabel, fontsize=25)
ax.set_ylabel(ylabel, fontsize=25)
ax.tick_params(axis='both', labelsize=20)
return ax
def Plot2D(img, title='', xlabel='', ylabel='', **kwargs):
ax = plt.subplot(111)
img = ax.imshow(img, origin='lower', **kwargs)
ax.set_title(title, fontsize=30, y=1.02)
ax.set_xlabel(xlabel, fontsize=25)
ax.set_ylabel(ylabel, fontsize=25)
ax.tick_params(axis='both', labelsize=20)
return ax, img
def UglyLoader(filename):
hdul = fits.HDUList(fits.open(filename))
nfluxes = hdul[0].header['NFLUXES']
......@@ -321,6 +533,36 @@ def UglyLoader(filename):
if __name__ == '__main__':
fname = ('/home/peter/Dokumente/Uni/Paris/Stage/'
'FirstSteps/Completness/montecarlo_results/nocources/'
'180419_105221_nosources_thresh2.5_nsim1000.fits')
rsfname = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'sourcemap/180418_141500_sourcemap_thresh2.fits')
rsfname = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'sourcemap/180419_183430_sourcemap_thresh2.fits')
hdul = fits.open(fname)
sources = Table.read(hdul['DETECTED_SOURCES'])
realsources = Table.read(rsfname, 'DETECTED_SOURCES')
mask = hdul['MATCHMASK'].data != 0
wcs = WCS(hdul[0].header)
nsim = hdul[0].header['NSIM']
aa = PCAreaEvaluation(mask=mask, sources=sources, wcs=wcs, nsim=nsim,
realsources=realsources, threshold_range=(2.5, 10),
threshold_bins=50)
if(1): # plot mask
fig = plt.figure()
ax, img = aa.PlotMask(sources=True, realsources=True)
plt.show(block=True)
if(1):
plt.figure()
aa.PlotDarkcountsThresholds()
plt.show(block=True)
'''
DATA_DIR = "/home/peter/Dokumente/Uni/Paris/Stage/data/v_1"
data = NikaMap.read(Path(DATA_DIR) / '..' / 'map.fits')
......@@ -332,6 +574,8 @@ if __name__ == '__main__':
fname = ('/home/peter/Dokumente/Uni/Paris/Stage/'
'FirstSteps/Completness/NEW_combine_fct_result.fits')
'''
'''
FLUX, SOURCE, FSOURCE = UglyLoader(fname)
xx = PCEvaluation(SOURCE, FSOURCE, sh, wcs, FLUX, mapbins=19,
......@@ -352,13 +596,34 @@ if __name__ == '__main__':
ylabel='Purity')
plt.show(block=True)
'''
'''
fname = ('/home/peter/Dokumente/Uni/Paris/Stage/'
'FirstSteps/Completness/montecarlo_results/nocources/'
'180417_174539_nosources_thresh2.5_nsim5000.fits')
rsfname = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'sourcemap/180418_120641_sourcemap_thresh2.fits')
rsfname = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'sourcemap/180418_141500_sourcemap_thresh2.fits')
sources = Table.read(fname, 1)
zz = PCEvaluation(sources, None, sh, wcs, None, mapbins=19,
threshold_bins=6, threshold_range=(2.5, 5))
zz.PlotAbsoluteDetection(thresh=2.5, normalize=5000.)
realsources = Table.read(rsfname, 1)
zz = PCEvaluation(sources, None, shape=sh, wcs=wcs, flux=None, mapbins=2,
threshold_bins=7, threshold_range=(2.5, 7),
realsources=realsources)
zz.PlotAbsoluteDetection(thresh=3.5, normalize=5000.)
plt.figure()
zz.PlotDarkcountsThresholds(1, 1, normalize=5000)
plt.figure()
#plt.imshow(zz.real_sources_map[:, :, 0], origin='lower')
#print(zz.real_sources_map.shape)
plt.show(block=True)
'''
#print(sources.keys())
#print(mask)
# PCAreaEvaluation(sources, mask=np.)
......@@ -158,15 +158,22 @@ def purity_worker(shape, wcs, sources, max_threshold=2):
corresponding 2D :class:`numpy.ndarray`
"""
sources_snr = sources['SNR']
# As we are interested by the cumulative numbers, keep all inside the
# upper pixel
sources_snr[sources_snr > max_threshold] = max_threshold
xx, yy, zz = wcs.wcs_world2pix(sources['ra'], sources['dec'],
sources_snr, 0)
if sources is not None:
sources_snr = sources['SNR']
# As we are interested by the cumulative numbers, keep all inside the
# upper pixel
sources_snr[sources_snr > max_threshold] = max_threshold
xx, yy, zz = wcs.wcs_world2pix(sources['ra'], sources['dec'],
sources_snr, 0)
'''
print(zz.shape)
plt.plot(zz)
plt.show(block=True)
sys.exit()
'''
# Number of fake sources recovered
if 'fake_sources' in sources.keys():
if sources is not None and 'fake_sources' in sources.keys():
_purity, _ = np.histogramdd(np.asarray([xx, yy, zz]).T + 0.5,
bins=np.asarray(shape),
range=list(zip([0]*len(shape), shape)),
......@@ -176,11 +183,15 @@ def purity_worker(shape, wcs, sources, max_threshold=2):
_purity = np.cumsum(_purity[..., ::-1], axis=2)[..., ::-1]
else:
_purity = None
# Number of total detected sources at a given threshold
_norm_pur, _ = np.histogramdd(np.asarray([xx, yy, zz]).T + 0.5,
bins=np.asarray(shape),
range=list(zip([0]*len(shape), shape)))
_norm_pur = np.cumsum(_norm_pur[..., ::-1], axis=2)[..., ::-1]
if sources is not None:
# Number of total detected sources at a given threshold
_norm_pur, _ = np.histogramdd(np.asarray([xx, yy, zz]).T + 0.5,
bins=np.asarray(shape),
range=list(zip([0]*len(shape), shape)))
_norm_pur = np.cumsum(_norm_pur[..., ::-1], axis=2)[..., ::-1]
else:
_norm_pur = None
return _purity, _norm_pur
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment