Commit 510fb286 authored by LUSTIG Peter's avatar LUSTIG Peter

Added pure monte carlo detection function

parent 221a157d
from __future__ import absolute_import, division, print_function
from pathlib import Path
import os
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import Pool, cpu_count
from functools import partial
from astropy import units as u
from astropy.io import ascii
from astropy.wcs import WCS
from astropy.utils.console import ProgressBar
from astropy.table import vstack, Table
from scipy.optimize import curve_fit
from nikamap import NikaMap, Jackknife
from nikamap.utils import pos_uniform
from astropy.io import fits
import argparse
'''
%load_ext autoreload
%autoreload 2
%matplotlib tk
'''
plt.ion()
def fake_worker(jkiter, min_threshold=2, nsources=8**2, flux=1*u.Jy,
within=(0, 1), cat_gen=pos_uniform, **kwargs):
"""The completness purity worker, create a fake dataset from an jackknifed
image and return catalogs
Parameters
----------
img : :class:`nikamap.NikaMap`
Jackknifed dataset
min_threshold : float
minimum threshold for the detection
**kwargs :
arguments for source injection
"""
img = jkiter()
# img = IMAGE
# Renormalize the stddev
std = img.check_SNR()
img.uncertainty.array *= std
# Actually rather slow... maybe check the code ?
print(flux)
img.add_gaussian_sources(nsources=nsources, within=within, peak_flux=flux,
cat_gen=cat_gen, **kwargs)
# ... match filter it ...
mf_img = img.match_filter(img.beam)
std = mf_img.check_SNR()
# print(std)
mf_img.uncertainty.array *= std
# print(mf_img.wcs)
# ... and detect sources with the lowest threshold...
# The gaussian fit from subpixel=True is very slow here...
mf_img.detect_sources(threshold=min_threshold)
return mf_img.sources, mf_img.fake_sources
plt.close('all')
DATA_DIR_SERVER = "/data/NIKA/Reduced/G2_COMMON_MODE_ONE_BLOCK/v_1/"
DATA_DIR_MYPC = "/home/peter/Dokumente/Uni/Paris/Stage/data/v_1"
parser = argparse.ArgumentParser()
parser.add_argument("-r", "--raw",
action="store",
dest="rawdirectory", default="mypc",
help="Path of Raw File Directory or Key.")
parser.add_argument("-t", "--threshold",
action="store",
dest="threshold", type=float, default=3.,
help="SNR Detection Threshold")
parser.add_argument("-n", "--nsim", action="store",
dest="nsim", type=int, default=4,
help="Number of Simulations")
parser.add_argument("-c", "--cores", action="store",
dest="cores", type=int, default=1,
help="Number of Cores")
parser.add_argument("-f", "--flux", action="store",
dest="flux", type=float, default=10.,
help="Flux of the Created Sources in mJy")
parser.add_argument("-s", "--sources", action="store",
dest="nsources", type=int, default=64,
help="Number of Sources Injected per Simulation")
'''
parser.add_argument("-m", "--method", action="store",
dest="method", default="unknown",
help="Name of the Reduction Method for Outfile Header")
'''
parser.add_argument("-o", "--outfile", action="store",
dest="outfile", default=None,
help=("Output File Name. If not specified auto-generated"
"from parameters"))
parser.add_argument("-d", "--outdir", action="store",
dest="outdir", default="montecarlo_results/",
help="Output directory. Must exist.")
args = parser.parse_args()
raw = args.rawdirectory
nsim = args.nsim
ncores = args.cores
flux = args.flux * u.mJy
# method = args.method
outfile = args.outfile
outdir = args.outdir
nsources = args.nsources
min_detection_threshold = args.threshold
if raw == "mypc":
DATA_DIR = DATA_DIR_MYPC
elif raw == "server":
DATA_DIR = DATA_DIR_SERVER
else:
DATA_DIR = raw
if outfile is None:
outfile = ('flux{}mJy_thresh{}_nsim{}.fits'
.format(flux.to_value(unit=u.mJy),
min_detection_threshold,
nsim))
outfile = outdir + outfile
'''
IMAGE = NikaMap.read(('/home/peter/Dokumente/Uni/Paris/'
'Stage/N2CLS/ForMarseille/'
'HLS091828_common_mode_kids_out/map_JK.fits'))
IMAGE.plot()
plt.show(block=True)
print(IMAGE)
'''
message = '''\
###############################################################
# Running Monte Carlo Tests #
# #
# Number of Simulations: {:5d} #
# Number of CPUs used: {:5d} #
# Minimum Detection Threshold:{:5.1f} #
# Flux of Injected Sources: {:5.1f} #
# Number of Injected Sources: {:5d} #
###############################################################
'''.format(nsim, ncores, min_detection_threshold, flux, nsources)
print(message)
jk_filenames = list(Path(DATA_DIR).glob('*/map.fits'))
for i in range(len(jk_filenames)):
jk_filenames[i] = str(jk_filenames[i])
FAKE_SOURCES = []
DETECTED_SOURCES = []
jk_iter = Jackknife(jk_filenames, n=nsim)
jk_iter_list = [jk_iter] * nsim
# print(type(jk_iter()))
p = Pool(ncores)
helpfunc = partial(fake_worker, **{'min_threshold': min_detection_threshold,
'nsources': nsources, 'flux': flux,
'within': (0, 1), 'cat_gen': pos_uniform})
# sources, fake_sources = helpfunc(jk_iter_list)
print('Begin Monte Carlo')
if(1):
res = p.map(helpfunc, jk_iter_list)
res = list(zip(*res))
DETECTED_SOURCES, FAKE_SOURCES = res[:]
else:
DETECTED_SOURCES = []
FAKE_SOURCES = []
with ProgressBar(nsim) as bar:
for iloop in range(nsim):
tmpsource, tmpfakesource = helpfunc(jk_iter)
DETECTED_SOURCES.append(tmpsource)
FAKE_SOURCES.append(tmpfakesource)
bar.update()
# To merge all the fake_sources and sources catalogs
fake_sources = Table()
sources = Table()
for _fake, _detected in zip(FAKE_SOURCES[:], DETECTED_SOURCES[:]):
n_fake = len(fake_sources)
n_detected = len(sources)
if _detected is not None:
_detected['ID'] = _detected['ID'] + n_detected
_detected['fake_sources'] = _detected['fake_sources'] + n_fake
sources = vstack([sources, _detected])
_fake['ID'] = _fake['ID'] + n_fake
_fake['find_peak'] = _fake['find_peak'] + n_detected
fake_sources = vstack([fake_sources, _fake])
phdu = fits.PrimaryHDU()
phdu.header['influx'] = '{}'.format(flux)
phdu.header['nsim'] = nsim
phdu.header['sourcespersim'] = nsources
phdu.header['dthresh'] = min_detection_threshold
hdul = [phdu]
if len(sources) > 0:
hdul.append(fits.BinTableHDU(data=sources, name='Detected_Sources'))
hdul.append(fits.BinTableHDU(data=fake_sources, name='Fake_Sources'))
hdul = fits.HDUList(hdul)
hdul.writeto(outfile, overwrite=True)
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