Commit 02d0f213 authored by LUSTIG Peter's avatar LUSTIG Peter
Browse files

added files

parent 2397d998
import sys
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
from completness2d import BasicEvaluation, LevMarLSQFitter, ChiSquaredPerNDF
from nikamap import NikaMap
from astropy.io import fits
from astropy.table import Table, vstack
import astropy.units as u
import script_purity as purity
import script_fluxboost as completeness
import numpy as np
from fitmodels import FluxBoosting
import matplotlib.pyplot as plt
import pynverse as pnvs
from pathlib import Path
import sys
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/data/test')
from picklesave import LoadObject
from counts_2mm import AbsoluteCountsCorrector, CombineCounts
def SaveFigure(figure, filename, outdir='/home/peter/plots'):
outpath = Path(outdir)
formats = ['.png', '.svg', '.pdf']
for format in formats:
_filename = outpath / (filename + format)
figure.savefig(str(_filename))
'''
class AbsoluteCountsCorrector(BasicEvaluation):
def __init__(self, realcounts, completeness, falsecounts):
assert np.array_equal(realcounts.mask, completeness.mask), (
"realcounts do not "
"cover the same area "
"as completeness")
assert np.array_equal(realcounts.mask, falsecounts.mask), (
"realcounts do not "
"cover the same area "
"as falsecounts")
self.realsources = realcounts.sources[~realcounts.sourceslocmask]
self.completeness = completeness
self.falsesources = falsecounts.sources[~falsecounts.sourceslocmask]
self.nfalsesim = falsecounts.no_sim
self.mask = realcounts.mask
self.shape = realcounts.shape
self.wcs = realcounts.wcs
(self.completenessfct, _), _ = completeness.FitCompleteness(0, weighted=True)
self.Deboost()
self.corrected_counts, _,_ ,_ = self.Count()
def Deboost(self):
completeness = self.completeness
realsources = self.realsources
falsesources = self.falsesources
(boostfct, _), _ = completeness.FitFluxboost(relative=False)
invboost = pnvs.inversefunc(boostfct)
realsources['deboosted_flux'] = invboost(
u.Quantity(realsources['flux_psf']).to(u.mJy))
falsesources['deboosted_flux'] = invboost(
u.Quantity(falsesources['flux_psf']).to(u.mJy))
return
def Count(self, bins=[0, 10]):
realsourcesflux = self.realsources['deboosted_flux'].copy()
falsesourcesflux = self.falsesources['deboosted_flux'].copy()
completenessfct = self.completenessfct
nfalsenim = self.nfalsesim
cp = self.completeness
realsourcesflux[realsourcesflux > bins[-1]] = bins[-1]
falsesourcesflux[falsesourcesflux > bins[-1]] = bins[-1]
realhist, edges = np.histogram(realsourcesflux, bins=bins)
# realhist = np.cumsum(realhist[::-1])[::-1]
falsehist= np.histogram(falsesourcesflux, bins=bins)[0] / nfalsenim
# falsehist = np.cumsum(falsehist[::-1])[::-1]
centers = (edges[0:-1] + edges[1:]) / 2
norm = completenessfct(centers)
norm = completenessfct(edges[0:-1])
corrected_counts = (realhist - falsehist) / norm
return corrected_counts, realhist, falsehist, norm
def CombineCounts(corrected_counts, completenesses, surfaces,
completeness_min=.5):
nbins = len(corrected_counts[0])
surfaces = np.tile(surfaces, (nbins, 1)).T
cpmask = completenesses < completeness_min
corrected_counts = np.ma.masked_array(corrected_counts, mask=cpmask)
surfaces = np.ma.masked_array(surfaces, mask=cpmask)
ccounts = corrected_counts.sum(axis=0)
addsurf = surfaces.sum(axis=0)
comb = corrected_counts.sum(axis=0) / surfaces.sum(axis=0)
return comb, ccounts, addsurf, cpmask, addsurf
'''
allmasks = purity.snrmasks
hlsmask = purity.hlsmask
bins = np.linspace(0, 10, 10)
bins = np.linspace(0, 10, 20)
# bins = np.geomspace(0.01, 10, 10)
bincenters = (bins[1:] + bins[0:-1]) / 2
ac = AbsoluteCountsCorrector(purity.rc, completeness.cp_phot, purity.fc)
completeness.cp_phot.FitFluxboost()
corr_counts, real_counts, false_counts, compl = ac.Count(bins=bins)
newtable = ac.realsources
#corr_counts * 60**2
corr_counts
real_counts
false_counts
compl
# counts = [ac.Count(bins=bins)]
surfaces = [ac.Surface()]
for i in range(1, 8):
print('mask', i)
_mask = allmasks[i] | hlsmask
purity.rc.ChangeMask(_mask)
completeness.cp_phot.ChangeMask(_mask)
purity.fc.ChangeMask(_mask)
ac = AbsoluteCountsCorrector(purity.rc, completeness.cp_phot, purity.fc)
# counts.append(ac.Count())
t1, t2, t3, t4 = ac.Count(bins=bins)
corr_counts = np.vstack((corr_counts, t1))
real_counts = np.vstack((real_counts, t2))
false_counts = np.vstack((false_counts, t3))
compl = np.vstack((compl, t4))
newtable = vstack([newtable, ac.realsources])
surfaces.append(ac.Surface())
# np.sum(newtable['completeness'] > .8)
# np.amax(newtable['flux_psf'])
here = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
# newtable.write(here / 'Catalog1mm.fits', overwrite=True)
print('Done')
plt.show()
len(surfaces)
surfaces = u.Quantity(surfaces)
sunit = u.arcmin**2
print(sunit)
surfaces = np.array(surfaces.to_value())
res, ccounts, surfs, cpmask, surfacesnorm = CombineCounts(corr_counts,
compl, surfaces,
completeness_min=.9)
# comb_err = ccounts /
res
res
# ccounts*60**2
resmask = res.mask
# res = res / sunit
print(res)
res = res.filled(np.nan) * (1/sunit)
type(res)
res = res.to(1/u.deg**2)
res
corr_counts.shape
compl
np.tile(surfaces, (9, 1)).T
real_counts
absdet = np.sum(np.ma.masked_array(real_counts, mask=cpmask), axis=0)
surfacesnorm
# print(sunit)
ccounts*60**2
np.sqrt(absdet)
ccounts / np.sqrt(absdet)
ccounts
errres = ccounts / np.sqrt(absdet) / surfacesnorm / sunit
errres
print(errres.to(1/u.deg**2))
errres.to(1/u.deg**2)
absdet = absdet / sunit
absdet
#%% Plot
%matplotlib tk
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/matthieu/Analyse/')
import number_counts_plots as ncp
ncp.ax2
import importlib
importlib.reload(ncp)
axdiff = ncp.ax2
ax=ncp.ax
# %%
tmp = np.cumsum(res[::-1])[::-1].value
tmperr = errres.to_value(1/u.deg**2)
tmperr = np.sqrt(np.cumsum(tmperr[::-1]**2)[::-1])
tmperr
tmp
absdet
tmperr
plotidx = tmp > 0
oldleg = ax.legend(fontsize=15, loc='lower left')
#ax.scatter(bincenters, tmp, color='r')
#ax.scatter(bins[:-1], tmp, color='y')
newhandle = ax.errorbar(bins[:-1][plotidx], tmp[plotidx], tmperr[plotidx],
label='New', color='r', linestyle='None', marker='x')
ax.set_xlabel('F$_{1200}$ [mJy]', fontsize=20)
ax.set_ylabel('N($>$F$_{1200}$) [deg$^{-2}$]', fontsize=20)
ax.tick_params('both', labelsize=15)
ax.set_title('Integral Counts', y=1.02, fontsize=25)
# newhandle
if 1:
CosmicVarianceDir = Path('/home/peter/Dokumente/Uni/Paris/Stage/CosmicVariance')
sidesvariance = LoadObject(str(CosmicVarianceDir / 'counts_SNIKA1200_new.pkl'))
intcounts = np.cumsum(res[::-1])[::-1]
tmperr = errres.to(1/u.deg**2)
tmperr = np.sqrt(np.cumsum(tmperr[::-1]**2)[::-1])
tmperr
ax = sidesvariance.Plot(ax=ax)
type(intcounts)
type(tmperr)
type(bins)
#ax.errorbar(bins[:-1], intcounts.to_value(1/u.deg**2)[:],
# tmperr.to_value(1/u.deg**2)[:], linestyle='None',
# marker='x', color='r', label='HLS Field')
ax.legend(fontsize=20)
ax.set_xlabel(r'$F$ [mJy]', fontsize=20)
ax.set_title('Integral Counts 1mm', fontsize=25, y=1.02)
# ax.set_xlim(xmax=3)
if 0:
import integral1mm
sidesbins = integral1mm.bins[:-1]
sidescounts = integral1mm.Nint
ax.scatter(sidesbins, sidescounts)
newhandle
leghandle0 = ax.plot([np.nan, np.nan], [np.nan, np.nan], color='blue')
leghandle1 = ax.fill(np.nan, np.nan, 'orange')
leghandle2 = ax.fill(np.nan, np.nan, 'yellow')
ax.legend([newhandle, (leghandle2[0], leghandle1[0], leghandle0[0])], ['New', 'SIDES'], loc='upper right', fontsize=15)
ax.add_artist(oldleg)
fig = ax.get_figure()
fig.set_size_inches(7, 5)
fig.subplots_adjust(left=0.15, bottom=0.15, right=0.98)
ax.set_xlim(xmax=12)
# SaveFigure(fig, '1mm_counts')
#%%
# %matplotlib tk
plotidx = absdet.to_value() > 1
len(plotidx)
len(res)
len(errres)
axdiff = ncp.ax2
print(res)
xsides = (sidesvariance.bins[:-1] + sidesvariance.bins[1:]) / 2
sidesds = (sidesvariance.bins[1:] - sidesvariance.bins[:-1])
myds = bins[1:] - bins[:-1]
myds
absdet
# sidesvariance.radius
# ysides = np.mean(sidesvariance.counts, axis=0) / (sidesvariance.radius**2 * np.pi)
ysides = sidesvariance.totalcounts
yerrsides = np.std(sidesvariance.counts, axis=0) / (sidesvariance.radius**2 * np.pi)
sidesvariance.counts.shape
oldleg = axdiff.legend(fontsize=15, loc='lower left')
# sidesvariance.totalcounts
# axdiff.scatter(bincenters, res.to_value(1/u.sr) * (bincenters*1E-3)**2.5*1000, color='yellow')
newhandle = axdiff.errorbar(bincenters[plotidx], (res.to_value(1/u.sr) * (bincenters*1E-3)**2.5*1000/myds)[plotidx],
(errres.to_value(1/u.sr) * (bincenters*1E-3)**2.5*1000/myds)[plotidx], color='red',
linestyle='None', marker='x')
# axdiff.plot(xsides, ysides.to_value(1/u.sr) * (xsides.to_value(u.Jy))**2.5*1000)
ysides = ysides.to_value(1/u.sr) * (xsides.to_value(u.Jy))**2.5 / sidesds.to_value(u.Jy)
yerrsides = yerrsides.to_value(1/u.sr) * (xsides.to_value(u.Jy))**2.5 / sidesds.to_value(u.Jy)
axdiff.plot(xsides, ysides, color='k')
axdiff.fill_between(xsides, ysides-2*yerrsides, ysides+2*yerrsides, color='yellow')
axdiff.fill_between(xsides, ysides-yerrsides, ysides+yerrsides, color='orange')
leghandle0 = ax.plot([np.nan, np.nan], [np.nan, np.nan], color='blue')
leghandle1 = ax.fill(np.nan, np.nan, 'orange')
leghandle2 = ax.fill(np.nan, np.nan, 'yellow')
axdiff.legend([newhandle, (leghandle2[0], leghandle1[0], leghandle0[0])], ['New', 'SIDES'], loc='upper right', fontsize=15)
axdiff.add_artist(oldleg)
axdiff.plot(xsides, ysides)
axdiff.set_xlabel(r'$F$ [mJy]', fontsize=20)
axdiff.set_ylabel(r'$\mathrm{d}N/\mathrm{d}F\cdot F^{5/2}$ [Jy$^{1.5}$/sr]', fontsize=20)
figdiff = axdiff.get_figure()
figdiff.subplots_adjust(bottom=0.15, left=.15, right=0.97)
ylim = axdiff.get_ylim()
axdiff.set_ylim(ylim)
axdiff.tick_params('both', labelsize=15)
#SaveFigure(figdiff, 'DifferentialCounts1mmWithSides')
#axdiff.plot(xsides, ysides.to_value(1/u.sr) * (xsides.to_value(u.Jy))**2.5 / sidesds.to_value(u.Jy))
plt.show()
print(res*bincenters**2.5)
# %%
tmp = counts[0] / ac.Surface()
tmp
tmp.to(1/u.deg**2)
np.array(tmp)
tmp = np.cumsum(tmp[::-1])[::-1]
len(bincenters)
len(tmp)
# %matplotlib tk
plt.plot(bincenters, tmp.to_value(1 / u.deg**2))
plt.gca().set_yscale('log')
# %% Sidescounts
CosmicVarianceDir = Path('/home/peter/Dokumente/Uni/Paris/Stage/CosmicVariance')
sidesvariance = LoadObject(str(CosmicVarianceDir / 'counts_SNIKA1200.pkl'))
ax = sidesvariance.Plot()
#Scott K S at al. (2012, AzTEC combined)
corr1100to1200 = 1. / 1.25 #(from SED at z=2)
S1100_s12 = np.array([1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.])*corr1100to1200
NsupS1100_s12 = np.array([1890.,750.,330.,150.,67.,32.,17.,9.,5.5,3.3,2.1])
NsupS1100l_s12 = np.array([120.,50.,30.,20.,13.,8.,4.,3.,2.1,1.6,1.8])
NsupS1100u_s12 = np.array([110.,50.,30.,20.,12.,6.,4.,2.,1.7,1.4,1.1])
ax.errorbar(S1100_s12, NsupS1100_s12,
yerr = (NsupS1100l_s12,NsupS1100u_s12),
fmt = 'g<',
label = 'Scott et al. (SD)')
ax.errorbar(bins[:-1][plotidx], tmp[plotidx], tmperr[plotidx],
label='HLS Field', color='r', linestyle='None', marker='x')
ax.set_xlabel(r'$F$ [mJy]', fontsize=20)
ax.set_title('Integral Counts 1mm', fontsize=25, y=1.02)
ax.set_xlim(xmax=13)
ax.legend(fontsize=20)
# %%
# SaveFigure(ax.get_figure(), 'IntegralCounts1mmWithSides')
newtable
if __name__ == '__main__' and 1:
diffcounts = res
plotmask = diffcounts > 0
differr = errres.to(1/u.deg**2)
ax = sidesvariance.Plot(integral=False)
ax.errorbar((bins[:-1])[plotmask], diffcounts.to_value(1/u.deg**2)[plotmask],
differr.to_value(1/u.deg**2)[plotmask], linestyle='None',
marker='x', color='r', label='HLS Field')
ax.legend(fontsize=20)
ax.set_xlabel(r'$F$ [mJy]', fontsize=20)
ax.set_title('Differential Counts 1mm', fontsize=25, y=1.02)
ax.set_ylabel(r'$\mathrm{d}N/\mathrm{d}F$', fontsize=20)
ax.set_xlim(xmax=3)
# SaveFigure(ax.get_figure(), 'DifferentialCounts1mmWithSides')
import sys
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
from completness2d import BasicEvaluation, LevMarLSQFitter, ChiSquaredPerNDF
from nikamap import NikaMap
from astropy.io import fits
from astropy.table import Table, vstack
import astropy.units as u
import numpy as np
from fitmodels import FluxBoosting
import matplotlib.pyplot as plt
import pynverse as pnvs
from pathlib import Path
# import integral2mm
import importlib
import sys
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/data/test')
from picklesave import LoadObject
'''
nm = NikaMap.read('/home/peter/Dokumente/Uni/Paris/Stage/data/map.fits')
be = BasicEvaluation(Table(), nm.shape, nm.wcs, nm.mask)
be.Surface()
nm.surface().to(u.arcminute**2)
'''
def SaveFigure(figure, filename, outdir='/home/peter/plots'):
outpath = Path(outdir)
formats = ['.png', '.svg', '.pdf']
for format in formats:
_filename = outpath / (filename + format)
figure.savefig(str(_filename))
class CatalogCreator(BasicEvaluation):
def __init__(self, realcounts, completeness, falsecounts):
assert np.array_equal(realcounts.mask, completeness.mask), (
"realcounts do not "
"cover the same area "
"as completeness")
assert np.array_equal(realcounts.mask, falsecounts.mask), (
"realcounts do not "
"cover the same area "
"as falsecounts")
self.realsources = realcounts.sources[~realcounts.sourceslocmask]
self.completeness = completeness
self.falsesources = falsecounts.sources[~falsecounts.sourceslocmask]
self.nfalsesim = falsecounts.no_sim
self.mask = realcounts.mask
self.shape = realcounts.shape
self.wcs = realcounts.wcs
(self.completenessfct, _), _ = completeness.FitCompleteness(
0, weighted=True)
self.realsources, self.falsesources = self.Deboost()
self.corrected_counts, _, _, _ = self.Count()
def Deboost(self):
completeness = self.completeness
realsources = self.realsources
falsesources = self.falsesources
(boostfct, _), _ = completeness.FitFluxboost(relative=False)
invboost = pnvs.inversefunc(boostfct)
realsources['deboosted_flux'] = invboost(
u.Quantity(realsources['flux_psf']).to(u.mJy))
falsesources['deboosted_flux'] = invboost(
u.Quantity(falsesources['flux_psf']).to(u.mJy))
return realsources, falsesources
def Derive(callable, x, h=1E-5):
return (callable(x+h)-callable(x-h)) / (2*h)
class AbsoluteCountsCorrector(BasicEvaluation):
def __init__(self, realcounts, completeness, falsecounts):
assert np.array_equal(realcounts.mask, completeness.mask), (
"realcounts do not "
"cover the same area "
"as completeness")
assert np.array_equal(realcounts.mask, falsecounts.mask), (
"realcounts do not "
"cover the same area "
"as falsecounts")
self.realsources = realcounts.sources[~realcounts.sourceslocmask]
self.completeness = completeness
self.falsesources = falsecounts.sources[~falsecounts.sourceslocmask]
self.nfalsesim = falsecounts.no_sim
self.mask = realcounts.mask
self.shape = realcounts.shape
self.wcs = realcounts.wcs
(self.completenessfct, _), _ = completeness.FitCompleteness(0, weighted=True)
self.AddSourceCompleteness()
self.Deboost()
self.corrected_counts, _, _, _ = self.Count()
def Deboost(self):
completeness = self.completeness
realsources = self.realsources
falsesources = self.falsesources
(boostfct, _), _ = completeness.FitFluxboost(relative=False)
invboost = pnvs.inversefunc(boostfct)
realsources['deboosted_flux'] = invboost(
u.Quantity(realsources['flux_psf']).to(u.mJy))
realsources['deboosted_flux'].unit = u.mJy
falsesources['deboosted_flux'] = invboost(
u.Quantity(falsesources['flux_psf']).to(u.mJy))
falsesources['deboosted_flux'].unit = u.mJy
sim_unc = completeness.flux_unc
sim_flux = completeness.flux_centers.to_value(u.mJy)
# sim_unc = completeness.fluxmean_absolute_unc
sim_unc = completeness.flux_absolute_unc
realsources['flux_unc'] = np.interp(realsources['flux_psf'],
sim_flux, sim_unc)
realsources['flux_unc'].unit = u.mJy
propagated_unc = Derive(invboost, realsources['flux_psf'])
# print('', propagated_unc)
propagated_unc = np.abs(propagated_unc * realsources['eflux_psf'])
realsources['propagated_unc'] = propagated_unc
realsources['propagated_unc'].unit = u.mJy
def Count(self, bins=[0, 10]):
realsourcesflux = self.realsources['deboosted_flux'].copy()
falsesourcesflux = self.falsesources['deboosted_flux'].copy()
completenessfct = self.completenessfct
nfalsenim = self.nfalsesim
cp = self.completeness
realsourcesflux[realsourcesflux > bins[-1]] = bins[-1]
falsesourcesflux[falsesourcesflux > bins[-1]] = bins[-1]
realhist, edges = np.histogram(realsourcesflux, bins=bins)
# realhist = np.cumsum(realhist[::-1])[::-1]
falsehist = np.histogram(falsesourcesflux, bins=bins)[0] / nfalsenim
# falsehist = np.cumsum(falsehist[::-1])[::-1]
centers = (edges[0:-1] + edges[1:]) / 2
norm = completenessfct(centers)
norm = completenessfct(edges[0:-1])
corrected_counts = (realhist - falsehist) / norm
return corrected_counts, realhist, falsehist, norm
def AddSourceCompleteness(self):
realsources = self.realsources
completenessfct = self.completenessfct
completenessval = completenessfct(realsources['flux_psf'])
realsources['completeness'] = completenessval
def CombineCounts(corrected_counts, completenesses, surfaces,
completeness_min=.5):
nbins = len(corrected_counts[0])
surfaces = np.tile(surfaces, (nbins, 1)).T
cpmask = completenesses < completeness_min
corrected_counts = np.ma.masked_array(corrected_counts, mask=cpmask)
surfaces = np.ma.masked_array(surfaces, mask=cpmask)
ccounts = corrected_counts.sum(axis=0)
addsurf = surfaces.sum(axis=0)
comb = corrected_counts.sum(axis=0) / surfaces.sum(axis=0)
return comb, ccounts, addsurf, cpmask, addsurf
if __name__ == '__main__':
import script_purity_2mm as purity
import script_fluxboost_2mm as completeness
allmasks = purity.snrmasks
hlsmask = purity.hlsmask
bins = np.geomspace(0.01, 1.3, 20)
bins=np.linspace(0, 1.3, 5)
bincenters = (bins[1:] + bins[0:-1]) / 2
ac = AbsoluteCountsCorrector(purity.rc, completeness.cp_phot, purity.fc)
corr_counts, real_counts, false_counts, compl = ac.Count(bins=bins)
surfaces = [ac.Surface()]
purity.rc.PlotMask(purity.rc.sources[~purity.rc.sourceslocmask])
print(ac.realsources)
print(ac.realsources['completeness'])
newtable = ac.realsources
for i in range(1, 6):
print('mask', i)
_mask = allmasks[i] | hlsmask
purity.rc.ChangeMask(_mask)
completeness.cp_phot.ChangeMask(_mask)
purity.fc.ChangeMask(_mask)
purity.rc.PlotMask(purity.rc.sources[~purity.rc.sourceslocmask])
ac = AbsoluteCountsCorrector(purity.rc, completeness.cp_phot, purity.fc)
# counts.append(ac.Count())
t1, t2, t3, t4 = ac.Count(bins=bins)
corr_counts = np.vstack((corr_counts, t1))
real_counts = np.vstack((real_counts, t2))
false_counts = np.vstack((false_counts, t3))
compl = np.vstack((compl, t4))
surfaces.append(ac.Surface())
# print(ac.realsources['deboosted_flux'], ac.realsources['flux_psf'])
print(ac.realsources)
print(ac.realsources['completeness'])
newtable = vstack([newtable, ac.realsources])
plt.show()
here = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
# newtable.write(here / 'Catalog2mm_small.fits', overwrite=True)
newtable
print('Done')