Commit e1c8ee84 authored by LUSTIG Peter's avatar LUSTIG Peter
Browse files

just saving...

parent 4297dc26
......@@ -5,8 +5,15 @@ import astropy.units as u
import matplotlib.pyplot as plt
from astropy.wcs import WCS
import warnings
import sys
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness')
from fitmodels import ModifiedErrorFunction
from astropy.modeling.fitting import LevMarLSQFitter
plt.close('all')
def GetLocationMask(sources, mask, addmask=None):
if addmask is None:
addmask = np.zeros(len(sources), dtype=bool)
......@@ -28,7 +35,7 @@ def GetInsideMask(sources, shape, within=(0, 1)):
nout = np.sum(outside)
if nout != 0:
warnings.warn("{} sources detected outside shape".format(nout),
UserWarning)
UserWarning)
return outside
......@@ -166,7 +173,11 @@ class CompletnessArea(BasicEvaluation):
completness, norm_completness = Completness2D(
sources, fake_sources, threshold_edges, flux_edges,
fake_sourceslocmask)
self.completness = completness / norm_completness[:, np.newaxis]
normed_completness = completness / norm_completness[:, np.newaxis]
self.completness = normed_completness
self.completness_err = np.sqrt(np.absolute(
((normed_completness * (1 - normed_completness))
/ norm_completness[:, np.newaxis])))
def InterpolateThresholds(self, value):
completness = self.completness
......@@ -191,6 +202,61 @@ class CompletnessArea(BasicEvaluation):
return constfluxinterp, constthreshinterp
def FitCompletness(self, idx, constant='thresh', fitmodel=None,
weighted=False):
if constant == 'thresh':
y = self.completness[:, idx]
y_err = self.completness_err[:, idx]
x = self.flux_centers.to_value(u.mJy)
elif constant == 'flux':
y = self.completness[idx]
y_err = self.completness_err[idx]
x = self.threshold_edges[:-1]
else:
raise UserWarning('you must chose \'flux\' or thresh for constant')
# define model and fitter
if fitmodel is None:
x0_guess = np.interp(.5, y, x)
fitmodel = ModifiedErrorFunction(x0_guess, 1, 1, .5)
fitmodel.amplitude.fix = True
fitmodel.offset.fix = True
fitter = LevMarLSQFitter()
# fit
if weighted:
res = fitter(fitmodel, x, y, weights=1./y_err)
else:
res = fitter(fitmodel, x, y)
# create plot
# Layout
fig = plt.figure(figsize=(7, 6))
left, right, bottom, top = 0.16, 0.97, 0.15, 0.9
fraction = 1./3
ax1 = fig.add_axes([left, bottom + (top-bottom) * fraction,
right-left, (1-fraction) * (top-bottom)])
ax2 = fig.add_axes([left, bottom, right-left,
(top-bottom) * fraction], sharex=ax1)
# plot
ax1 = self.PlotCompletness1D(idx, constant=constant, marker='o',
ax=ax1)
plotx = np.linspace(x[0], x[-1], 100)
ax1.plot(plotx, res(plotx))
ax2.scatter(x, y - res(x))
ax2.axhline(0)
# modify labels
ax1.tick_params(axis='both', labelsize=15)
ax2.tick_params(axis='both', labelsize=15)
ax1.set_ylabel('Comppletness', fontsize=20)
ax2.set_ylabel('Comp-Fit', fontsize=20)
return res, [ax1, ax2]
def PlotCompletness1D(self, idx, constant='flux', linestyles=None,
ax=None, **kwargs):
......@@ -200,12 +266,14 @@ class CompletnessArea(BasicEvaluation):
if constant == 'thresh':
completnesses = self.completness[:, idx].T
completnesses_err = self.completness_err[:, idx].T
x = self.flux_centers.to_value(u.mJy)
xlabel = 'Flux [mJy]'
labelval = self.threshold_edges[idx]
legendtitle = 'SNR'
elif constant == 'flux':
completnesses = self.completness[idx]
completnesses_err = self.completness_err[idx]
x = self.threshold_edges[:-1]
xlabel = 'Detection Threshold'
labelval = self.flux_centers[idx].to_value(u.mJy)
......@@ -221,12 +289,13 @@ class CompletnessArea(BasicEvaluation):
fig.subplots_adjust(left=0.16, right=0.97, bottom=0.15, top=0.9)
ax = fig.add_subplot(111)
for i, (completness, linestyle) in enumerate(
zip(completnesses, linestyles)):
ax.plot(x, completness, linestyle=linestyle,
label='{:.2f}'.format(labelval[i]), **kwargs)
for i, (completness, cerror, linestyle) in enumerate(
zip(completnesses, completnesses_err,
linestyles)):
ax.errorbar(x, completness, cerror, linestyle=linestyle,
label='{:.2f}'.format(labelval[i]), **kwargs)
legend = ax.legend(fontsize=22, title=legendtitle, loc='upper right',
legend = ax.legend(fontsize=22, title=legendtitle, loc='lower right',
framealpha=1)
plt.setp(legend.get_title(), fontsize=20)
ax.set_title('Completness', fontsize=30, y=1.02)
......@@ -338,8 +407,9 @@ class FluxArea(CompletnessArea):
super().__init__(sources, fake_sources, shape, wcs, threshold_edges,
flux_edges, flux_centers, mask)
self.fluxpercentiles, self.snrbins = self.FluxPercentiles()
def FluxPercentiles(self, flux, snrbins=None,
def FluxPercentiles(self, snrbins=None,
percentiles=np.array([2.3, 15.9, 50., 84.1, 97.7])):
fake_sources = self.fake_sources
sources = self.sources
......@@ -354,19 +424,53 @@ class FluxArea(CompletnessArea):
xval[xval > bins[-1]] = bins[-1]
binidx = np.digitize(xval, bins, right=True)
percentiles = np.zeros((len(bins)-1, len(percentiles)))
outflux = sources[fake_sources['find_peak']]['flux_psf']
influx = fake_sources['amplitude']
perc = np.zeros((len(bins)-1, len(percentiles)))
outflux = u.Quantity(sources[fake_sources['find_peak'].filled(0)]
['flux_psf']).to_value(u.mJy)
influx = u.Quantity(fake_sources['amplitude']).to_value(u.mJy)
# print('outflux', type(outflux), len(outflux), outflux.unit, outflux[:5])
# print('influx', type(influx), len(influx), influx.unit, influx[:5])
fluxdeviation = np.ma.array(
(outflux - influx) / influx,
mask=fake_sources['find_peak'].mask)
for i in range(len(bins)-1):
binmask = binidx != i
totmask = binmask | fake_sources['find_peak'].mask | locmask
percentiles[i] = np.nanpercentile(fluxdeviation[~totmask],
percentiles)
tmp = np.array(fluxdeviation[~totmask])
if len(tmp) > 0:
perc[i] = np.nanpercentile(tmp, percentiles)
return perc, snrbins
def PlotPercentiles(self, exact=True, ax=None):
percentiles = np.copy(self.fluxpercentiles).T
medidx = percentiles.shape[0] // 2
flux_centers = self.flux_centers
if ax is None:
fig = plt.figure(figsize=(7, 6))
ax = fig.add_subplot(111)
return percentiles
if exact:
plt.plot(flux_centers, percentiles[medidx], linewidth=3,
label='Median')
colorarr = ['orange', 'yellow']
for iperc in range(medidx, 0, -1):
plt.fill_between(flux_centers,
percentiles[medidx-iperc],
percentiles[medidx+iperc],
alpha=1, facecolor=colorarr[iperc-1],
label=r'${{{}}}\sigma$'.format(iperc))
else:
raise NotImplementedError
plt.axhline(1, color='k')
plt.title('Flux Resolution', y=1.02, fontsize=25)
plt.xlabel('Injected Flux [mJy]', fontsize=20)
plt.ylabel(r'F$_{\mathrm{out}}$ / F$_{\mathrm{in}}$', fontsize=20)
plt.legend(fontsize=20, loc='upper right')
ax.set_yscale('log')
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
def Plot1D(x, y, title='', xlabel='', ylabel='', ax=None, **kwargs):
......@@ -461,10 +565,12 @@ class Purity:
ylabel=ylabel, ax=ax)
return ax
# %matplotlib tk
showcompletness = True
showfalsecounts = False
showrealcounts = False
showcompletness = False
showfalsecounts = True
showrealcounts = True
showflux = False
maskfilename = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/masks.fits')
......@@ -489,7 +595,7 @@ if showcompletness:
fluxedges = LogBinEdges(fake_sources['amplitude']) * u.Jy
fluxcenters = u.Quantity(np.squeeze(np.sort(np.unique(
fake_sources['amplitude']))))
# %%
plt.close('all')
# %matplotlib tk
obj = CompletnessArea(sources, fake_sources, shape=(501, 501), wcs=wcs,
......@@ -503,8 +609,6 @@ if showcompletness:
obj.PlotCompletness1D(np.array([1, 25, 49]), constant='thresh')
obj.PlotCompletness1D(np.array([10, 25, 30]), constant='flux')
obj.InterpolateThresholds(0.8)
# %%
if showfalsecounts:
emptycountsfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/'
......@@ -535,3 +639,36 @@ if showfalsecounts:
pur = Purity(abs_counts, false_counts, threshbins[0:-1])
pur.Plot()
if showflux:
filename = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/PhotometryNewrealisation3.5.fits')
with fits.open(filename) as hdul:
sources = Table.read(hdul, 'DETECTED_SOURCES')
fake_sources = Table.read(hdul, 'FAKE_SOURCES')
fluxedges = LogBinEdges(fake_sources['amplitude']) * u.Jy
fluxcenters = u.Quantity(np.squeeze(np.sort(np.unique(
fake_sources['amplitude']))))
threshbins = np.arange(3.5, 7, .5)
obj = FluxArea(sources, fake_sources, shape=(501, 501), wcs=wcs,
threshold_edges=threshbins, flux_edges=fluxedges,
mask=totmask, flux_centers=fluxcenters)
obj.FluxPercentiles()
obj.Plot2D()
# %matplotlib tk
idx = 0
ax = obj.PlotCompletness1D(idx, constant='thresh')
f, axes = obj.FitCompletness(idx=idx, constant='thresh')
# plt.show()
'''
f = ModifiedErrorFunction(1, 1, 1, .5)
f.amplitude.fix = True
f.offset.fix = True
fittedfunc, axes = obj.FitCompletness(0, fitmodel=f)
'''
print('done.')
# %%
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