Commit c5968bb1 authored by LUSTIG Peter's avatar LUSTIG Peter

End of the week...

parent 70724fe5
......@@ -105,8 +105,11 @@ def Combine(directory, key, outfile, precision=0.001*u.mJy):
directory = 'montecarlo_results/700/'
outname = 'NEW_combine_fct_result.fits'
key = '*700.fits'
directory = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'montecarlo_results/photometry/thresh4')
outname = 'allresults_phot_thresh4.fits'
key = '*.fits'
Combine(directory, key, outname)
......
This diff is collapsed.
......@@ -25,6 +25,9 @@ from mpl_toolkits.axes_grid1 import make_axes_locatable
import dill as pickle
from matplotlib.ticker import FormatStrFormatter
from collections import OrderedDict
import warnings
from astropy.modeling.models import Gaussian1D
from astropy.modeling.fitting import LevMarLSQFitter
def find_nearest(array, values):
......@@ -409,122 +412,95 @@ def CombineMeasurements(sourceslist, fakesourceslist):
return sources, fake_sources
if __name__ == "__main__":
DATA_DIR = "/home/peter/Dokumente/Uni/Paris/Stage/data/v_1"
data = NikaMap.read(Path(DATA_DIR) / '..' / 'map.fits')
bins = 19
flux_bins = 2
flux_range = [0.1, 10]
threshold_bins = 5
threshold_range = [2, 5]
threshold_range = [3, 7.5]
# Does not really make sense... better define edges
fluxes = np.logspace(np.log10(flux_range[0]), np.log10(flux_range[1]),
flux_bins)*u.mJy
fluxes_edges = np.logspace(np.log10(flux_range[0]), np.log10(flux_range[1]),
flux_bins + 1)*u.mJy
# Does not really make sense... better define edges
threshold = np.linspace(threshold_range[0], threshold_range[1], threshold_bins)
threshold_edges = np.linspace(threshold_range[0], threshold_range[1],
threshold_bins+1)
shape_4D, wcs_4D = completness_purity_wcs(data.shape, data.wcs, bins=bins,
flux_range=flux_range,
flux_bins=flux_bins, flux_log=True,
threshold_range=threshold_range,
threshold_bins=threshold_bins)
# Testing the lower edges
wcs_threshold = wcs_4D.sub([3])
assert np.all(np.abs(wcs_threshold.all_pix2world(np.arange(threshold_bins+1)-0.5, 0) - threshold_edges) < 1e-15)
wcs_flux = wcs_4D.sub([4])
assert np.all(np.abs(wcs_flux.all_pix2world(np.arange(flux_bins+1)-0.5, 0) - fluxes_edges.value) < 1e-13)
def pos_in_mask(pos, mask=None, nsources=1, extra=None, retindex=False):
"""Check if pos is in mask, issue warning with less than nsources
STOLEN FROM NIKAMAP AND MODIFIED TO RETURN INDEX
Parameters
----------
pos : array_like (N, 2)
pixel indexes (y, x) to be checked in mask
mask : 2D boolean array_like
corresponding mask
nsources : int
the requested number of sources
# DEBUG :
# flux, nsources, within, wcs, shape, nsim, jk_filenames = 10*u.mJy, 8**2,
# (0, 1), wcs_4D.sub([1, 2, 3]), shape_4D[0:3], np.multiply(*shape_4D[0:2])
# * 100, filenames
Returns
-------
:class:`numpy.ndarray`
the pixel indexes within the mask
"""
pos = np.asarray(pos)
inside = np.ones(pos.shape[0], dtype=bool)
# This is a single run check for a single flux
if mask is not None:
pos_idx = np.floor(pos + 0.5).astype(int)
inside = ~mask[pos_idx[:, 0], pos_idx[:, 1]]
if retindex:
return ~inside
pos = pos[inside]
hdul = fits.HDUList(fits.open('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/combined_tables_long.fits'))
# sys.exit()
nfluxes = hdul[0].header['NFLUXES']
print('{} different fluxes found'.format(nfluxes))
if pos.shape[0] < nsources:
warnings.warn("Only {} positions".format(pos.shape[0]), UserWarning)
# Get fluxlist:
if extra is None:
return pos
else:
return pos, extra[inside]
indata = []
for isimu in range(nfluxes):
_FLUX = u.Quantity(hdul[0].header['flux{}'.format(isimu)])
def Flux1D(detectedflux, catflux, bins=10, fitter=LevMarLSQFitter):
labelsize = 25
ticksize = 20
textfs = 25
boxprops = dict(facecolor='white', edgecolor='black',
boxstyle='round, pad=.3', alpha=.5)
_SOURCES = Table.read(hdul['DETECTED_SOURCES{}'
.format(_FLUX)])
_FAKE_SOURCES = Table.read(hdul['FAKE_SOURCES{}'
.format(_FLUX)])
indata.append([_FLUX, _SOURCES, _FAKE_SOURCES])
fluxrel = (detectedflux / catflux).decompose() - 1.
'''
flux = []
for isimu in range(nfluxes):
flux.append(u.Quantity(hdul[0].header['flux{}'.format(isimu)]))
print(fluxrel)
print(fluxrel.shape)
print(type(fluxrel))
print('creating histogram')
'''
# helpfunc = partial(Evaluate, **{'hdul': hdul})
p = Pool(cpu_count())
res = p.map(Evaluate, indata)
res = list(zip(*res))
FLUX = np.array(res[0])
COMPLETNESS = np.array(res[1])
PURITY = np.array(res[2])
idxsort = np.argsort(FLUX)
FLUX = FLUX[idxsort]
COMPLETNESS = COMPLETNESS[idxsort]
PURITY = PURITY[idxsort]
midbin = int(bins/2)
# sys.exit()
# %% PlotFigure
PlotEvaluation(COMPLETNESS[:, midbin, midbin, :], title='Completness',
flux=np.array(FLUX), thresh=threshold, cmap='bone',
nfluxlabels=10, nthreshlabels=5, aspect='auto')
PlotEvaluation(PURITY[:, midbin, midbin, :], title='Purity',
flux=np.array(FLUX), thresh=threshold, cmap='bone',
nfluxlabels=10, nthreshlabels=5, aspect='auto')
plt.show(block=True)
# %% 2D plot
PlotFixedThreshold(np.array((3, 5, 7)), (midbin, midbin), COMPLETNESS,
threshold, np.array(FLUX), nfluxlabels=None,
hlines=[.5, .8, .9])
# %% Plot map
# Plot_CompPur(completness, purity, threshold, nsim=None, savename=None,
# flux=None):
plotidx = np.array([24, 30, 40, 49])
# plotthreshidx = np.array([0, 4, 8])
plotcomp = COMPLETNESS[plotidx]
plotpur = PURITY[plotidx]
# plotthresh = threshold[plotthreshidx]
plotflux = np.array(FLUX[plotidx])
for i in range(len(plotidx)):
_comp = plotcomp[i]
_pur = plotpur[i]
Plot_CompPur(_comp, _pur, threshold, flux=plotflux[i])
plt.show(block=True)
title = 'Flux Resolution'
try:
iter(catflux)
# if it works maybe add fluxrange to title
except TypeError:
title += ' {:.2f}'.format(catflux)
plt.figure()
hist, edges, patches = plt.hist(fluxrel, bins=bins, density=True)
plt.title(title, fontsize=30, y=1.02)# , loc='left')
plt.xlabel(r'$\frac{F_{\mathrm{det}}}{F_{\mathrm{in}}}-1$',
fontsize=labelsize)
plt.ylabel('Normalized Counts', fontsize=labelsize)
plt.xticks(fontsize=ticksize)
plt.yticks(fontsize=ticksize)
center = edges[:-1] + (edges[1:] - edges[:-1]) / 2
finit = Gaussian1D(mean=0)
fitf = fitter()
f = fitf(finit, center, hist)
x = np.linspace(edges[0], edges[-1], 200)
ax = plt.gca()
plt.plot(x, f(x), color='r', linewidth=3)
fittext = ('Entries: {}\n'.format(len(detectedflux)) +
r'$x_0={:.2f}$'.format(f.mean.value) + '\n' +
r'$\sigma={:4.2f}$'.format(f.stddev.value))
# ax.text(1-0.05, 0.93, fittext,
ax.text(.5, 1-0.93, fittext,
fontsize=textfs, bbox=boxprops,
horizontalalignment='center',
verticalalignment='bottom',
ma='right',
transform=ax.transAxes)
return f
def DeletePercentiles(array, minclip=0, maxclip=100):
print(type(array))
percentiles = np.percentile(array, (minclip, maxclip))
print(percentiles)
mask = np.array((array < percentiles[0]) | (array > percentiles[1]),
dtype=bool)
return array[~mask]
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