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

corrected several errors...

parent e1c8ee84
......@@ -7,7 +7,7 @@ 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 fitmodels import ModifiedErrorFunction# , FluxRelativeBoosting
from astropy.modeling.fitting import LevMarLSQFitter
......@@ -134,25 +134,19 @@ class CompletnessArea(BasicEvaluation):
shape=None, wcs=None, threshold_edges=None, flux_edges=None,
flux_centers=None, mask=None):
sources = AddXYCoordinatesTo(sources, wcs)
fake_sources = AddXYCoordinatesTo(fake_sources, wcs)
if 'xmean' not in sources.keys():
sources = AddXYCoordinatesTo(sources, wcs)
if 'xmean' not in fake_sources.keys():
fake_sources = AddXYCoordinatesTo(fake_sources, wcs)
self.sources = sources
self.fake_sources = fake_sources
self.shape = shape
self.ChangeMask(mask)
self.threshold_edges = threshold_edges
self.flux_edges = flux_edges
self.flux_centers = flux_centers
self.Completness()
def ChangeMask(self, mask):
fake_sourceslocmask = np.zeros(len(fake_sources), dtype=bool)
sourceslocmask = np.zeros(len(sources), dtype=bool)
if mask is not None:
shape = self.shape
fake_sourceslocmask = GetLocationMask(fake_sources, mask)
# should check if corresponding fake_source is in mask or not
# sourceslocmask = fakesourceslocmask[sources['fake_sources']] ?
......@@ -162,6 +156,17 @@ class CompletnessArea(BasicEvaluation):
self.mask = mask
self.sourceslocmask = sourceslocmask
self.fake_sourceslocmask = fake_sourceslocmask
self.threshold_edges = threshold_edges
self.flux_edges = flux_edges
self.flux_centers = flux_centers
self.Completness()
def ChangeMask(self, mask):
self.__init__(self.sources, self.fake_sources,
shape=self.shape, wcs=None,
threshold_edges=self.threshold_edges,
flux_edges=self.flux_edges,
flux_centers=self.flux_centers, mask=mask)
def Completness(self):
threshold_edges = self.threshold_edges
......@@ -254,6 +259,7 @@ class CompletnessArea(BasicEvaluation):
ax1.set_ylabel('Comppletness', fontsize=20)
ax2.set_ylabel('Comp-Fit', fontsize=20)
ax2.set_xscale('log')
return res, [ax1, ax2]
......@@ -327,7 +333,10 @@ class CompletnessArea(BasicEvaluation):
flux[fluxlabelidx], decimals=2),
separator=',')[1:-1].split(',')
ax.imshow(completness, origin='lower', aspect='auto')
im = ax.imshow(completness, origin='lower', aspect='auto')
fig = ax.get_figure()
cbar = fig.colorbar(im)
cbar.ax.tick_params(labelsize=15)
ax.set_xticks(threshlabelidx)
ax.set_xticklabels(threshlabel, {'fontsize': 15})
......@@ -343,40 +352,38 @@ class CompletnessArea(BasicEvaluation):
class CountsAreaEvaluation(BasicEvaluation):
def __init__(self, sources,
shape=None, wcs=None, threshold_edges=None,
no_sim=1,
mask=None):
sources = AddXYCoordinatesTo(sources, wcs)
def __init__(self, sources, shape=None, wcs=None, threshold_edges=None,
no_sim=1, mask=None):
if 'xmean' not in sources.keys():
sources = AddXYCoordinatesTo(sources, wcs)
self.sources = sources
self.shape = shape
self.mask = mask
self.sourceslocmask = self.ChangeMask(mask)
self.threshold_edges = threshold_edges
self.no_sim = no_sim
self.counts = self.Counts()
def ChangeMask(self, mask):
sources = self.sources
sourceslocmask = np.zeros(len(sources), dtype=bool)
if mask is not None:
shape = self.shape
# should check if corresponding fake_source is in mask or not
# sourceslocmask = fakesourceslocmask[sources['fake_sources']] ?
sourceslocmask = GetInsideMask(sources, shape)
sourceslocmask = GetLocationMask(sources, mask,
addmask=sourceslocmask)
self.mask = mask
self.sourceslocmask = sourceslocmask
return sourceslocmask
self.threshold_edges = threshold_edges
self.no_sim = no_sim
self.counts = self.Counts()
def ChangeMask(self, mask):
self.__init__(sources=self.sources, shape=self.shape, wcs=None,
threshold_edges=self.threshold_edges, no_sim=self.no_sim,
mask=mask)
def Counts(self):
false_counts = AbsoluteCounts(self.sources, self.threshold_edges,
self.sourceslocmask)
false_counts /= self.no_sim
self.counts = false_counts
return false_counts
def PlotCounts(self, ax=None, label=None, **kwargs):
......@@ -409,6 +416,16 @@ class FluxArea(CompletnessArea):
flux_edges, flux_centers, mask)
self.fluxpercentiles, self.snrbins = self.FluxPercentiles()
def ChangeMask(self, mask):
self.__init__(sources=self.sources,
fake_sources=self.fake_sources,
shape=self.shape, wcs=None,
threshold_edges=self.threshold_edges,
flux_edges=self.flux_edges,
flux_centers=self.flux_centers,
mask=mask)
def FluxPercentiles(self, snrbins=None,
percentiles=np.array([2.3, 15.9, 50., 84.1, 97.7])):
fake_sources = self.fake_sources
......@@ -419,43 +436,115 @@ class FluxArea(CompletnessArea):
xval = u.Quantity(fake_sources['amplitude']).to_value(u.mJy)
if snrbins is not None:
bins = snrbins
xval = np.ma.array(sources[fake_sources['find_peak']]['SNR'],
xval = np.ma.array(sources[fake_sources['find_peak']
.filled(0)]['SNR'],
mask=fake_sources['find_peak'].mask)
xval[xval > bins[-1]] = bins[-1]
# xval[xval > bins[-1]] = bins[-1]
binidx = np.digitize(xval, bins, right=True)
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])
fake_sources['fluxgroup'] = - np.ones(len(fake_sources), dtype=int)
fluxdeviation = np.ma.array(
(outflux - influx) / influx,
outflux / influx,
mask=fake_sources['find_peak'].mask)
for i in range(len(bins)-1):
binmask = binidx != i
fake_sources['fluxdev'] = fluxdeviation
fake_sources['fluxdevgroup'] = binidx
for i in range(0, len(bins)-1):
binmask = binidx != (i + 1)
totmask = binmask | fake_sources['find_peak'].mask | locmask
tmp = np.array(fluxdeviation[~totmask])
if len(tmp) > 0:
if len(tmp) > 5:
perc[i] = np.nanpercentile(tmp, percentiles)
self.snrbins = snrbins
self.fluxpercentiles = perc
return perc, snrbins
'''
def FitFluxboost(self, model=None, asfunctionof='flux'):
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)
ax2.set_xscale('log')
return res, [ax1, ax2]
'''
def PlotPercentiles(self, exact=True, ax=None):
percentiles = np.copy(self.fluxpercentiles).T
medidx = percentiles.shape[0] // 2
flux_centers = self.flux_centers
x = self.flux_centers.to_value(u.mJy)
xlabel = 'Injected Flux [mJy]'
if self.snrbins is not None:
x = self.snrbins[:-1]
xlabel = 'SNR'
if ax is None:
fig = plt.figure(figsize=(7, 6))
ax = fig.add_subplot(111)
if exact:
plt.plot(flux_centers, percentiles[medidx], linewidth=3,
plt.plot(x, percentiles[medidx], linewidth=3,
label='Median')
colorarr = ['orange', 'yellow']
for iperc in range(medidx, 0, -1):
plt.fill_between(flux_centers,
plt.fill_between(x,
percentiles[medidx-iperc],
percentiles[medidx+iperc],
alpha=1, facecolor=colorarr[iperc-1],
......@@ -465,7 +554,7 @@ class FluxArea(CompletnessArea):
plt.axhline(1, color='k')
plt.title('Flux Resolution', y=1.02, fontsize=25)
plt.xlabel('Injected Flux [mJy]', fontsize=20)
plt.xlabel(xlabel, fontsize=20)
plt.ylabel(r'F$_{\mathrm{out}}$ / F$_{\mathrm{in}}$', fontsize=20)
plt.legend(fontsize=20, loc='upper right')
ax.set_yscale('log')
......@@ -568,9 +657,9 @@ class Purity:
# %matplotlib tk
showcompletness = False
showfalsecounts = True
showfalsecounts = False
showrealcounts = True
showflux = False
showflux = True
maskfilename = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/masks.fits')
......@@ -580,6 +669,9 @@ with fits.open(maskfilename) as hdul:
hlsmask = hdul['hls3fwhm'].data.astype(bool)
wcs = WCS(hdul['hitmask'].header)
anothermaskfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/verify/'
'snrmasks.fits')
hitmask = fits.getdata(anothermaskfile,'SNRMASK0').astype(bool)
totmask = hlsmask | hitmask
shape = totmask.shape
threshbins = np.arange(2.5, 7, .5)
......@@ -609,6 +701,7 @@ 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/'
......@@ -656,13 +749,17 @@ if showflux:
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()
# obj.FluxPercentiles()
obj.Plot2D(nthreshlabels=6)
# %matplotlib tk
idx = 0
ax = obj.PlotCompletness1D(idx, constant='thresh')
f, axes = obj.FitCompletness(idx=idx, constant='thresh')
# plt.show()
obj.PlotPercentiles()
obj.ChangeMask(totmask)
snrbins = np.linspace(3.5, 7, 8)
obj.FluxPercentiles(snrbins=snrbins)
obj.PlotPercentiles()
'''
f = ModifiedErrorFunction(1, 1, 1, .5)
......
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