Commit fe06153c authored by LUSTIG Peter's avatar LUSTIG Peter

Saving map masks, bug when parallel

parent 3766f21f
...@@ -23,6 +23,7 @@ from time import clock ...@@ -23,6 +23,7 @@ from time import clock
import argparse import argparse
import sys import sys
import datetime import datetime
from astropy.io.fits import ImageHDU
''' '''
%load_ext autoreload %load_ext autoreload
%autoreload 2 %autoreload 2
...@@ -34,7 +35,7 @@ plt.ion() ...@@ -34,7 +35,7 @@ plt.ion()
def fake_worker(jkiter, min_threshold=2, nsources=8**2, flux=1*u.Jy, def fake_worker(jkiter, min_threshold=2, nsources=8**2, flux=1*u.Jy,
within=(0, 1), cat_gen=pos_uniform, parity_threshold=1, within=(0, 1), cat_gen=pos_uniform, parity_threshold=1,
**kwargs): retmask=True, **kwargs):
"""The completness purity worker, create a fake dataset from an jackknifed """The completness purity worker, create a fake dataset from an jackknifed
image and return catalogs image and return catalogs
...@@ -70,7 +71,10 @@ def fake_worker(jkiter, min_threshold=2, nsources=8**2, flux=1*u.Jy, ...@@ -70,7 +71,10 @@ def fake_worker(jkiter, min_threshold=2, nsources=8**2, flux=1*u.Jy,
# The gaussian fit from subpixel=True is very slow here... # The gaussian fit from subpixel=True is very slow here...
mf_img.detect_sources(threshold=min_threshold) mf_img.detect_sources(threshold=min_threshold)
return mf_img.sources, mf_img.fake_sources if retmask:
return mf_img.sources, mf_img.fake_sources, img.mask, mf_img.mask
else:
return mf_img.sources, mf_img.fake_sources
plt.close('all') plt.close('all')
...@@ -81,12 +85,15 @@ DATA_DIR_MYPC = Path("/home/peter/Dokumente/Uni/Paris/Stage/data/v_1") ...@@ -81,12 +85,15 @@ DATA_DIR_MYPC = Path("/home/peter/Dokumente/Uni/Paris/Stage/data/v_1")
flux = np.geomspace(1, 10, 3) * u.mJy flux = np.geomspace(1, 10, 3) * u.mJy
flux = [0] flux = [0]
nsim = 2 flux = np.linspace(1, 3, 2)*u.mJy
flux = np.array([10]) * u.mJy
nsim = 4
min_detection_threshold = 3 min_detection_threshold = 3
nsources = 2
nsources = 0 nsources = 0
outdir = Path("montecarlo_results/") outdir = Path("montecarlo_results/")
outdir = Path("testdir") outdir = Path("testdir")
ncores = 2 ncores = 4
timeprefix = '{date:%y%m%d_%H%M%S}_'.format(date=datetime.datetime.now()) timeprefix = '{date:%y%m%d_%H%M%S}_'.format(date=datetime.datetime.now())
...@@ -133,6 +140,7 @@ print('Done in {:.2f}s'.format(clock()-t0)) ...@@ -133,6 +140,7 @@ print('Done in {:.2f}s'.format(clock()-t0))
print('Begin Monte Carlo') print('Begin Monte Carlo')
jk_iter_list = [jk_iter] * nsim jk_iter_list = [jk_iter] * nsim
p = Pool(ncores) p = Pool(ncores)
globmask = None
if nsources != 0: if nsources != 0:
for _flux in flux: for _flux in flux:
...@@ -151,15 +159,25 @@ if nsources != 0: ...@@ -151,15 +159,25 @@ if nsources != 0:
res = p.map(helpfunc, jk_iter_list) res = p.map(helpfunc, jk_iter_list)
res = list(zip(*res)) res = list(zip(*res))
DETECTED_SOURCES, FAKE_SOURCES = res[:] (DETECTED_SOURCES,
FAKE_SOURCES,
immask,
matchmask) = res
else: else:
DETECTED_SOURCES = [] DETECTED_SOURCES = []
FAKE_SOURCES = [] FAKE_SOURCES = []
immask = []
matchmask = []
with ProgressBar(nsim) as bar: with ProgressBar(nsim) as bar:
for iloop in range(nsim): for iloop in range(nsim):
tmpsource, tmpfakesource = helpfunc(jk_iter) (tmpsource,
tmpfakesource,
_immask,
_matchmask) = helpfunc(jk_iter)
DETECTED_SOURCES.append(tmpsource) DETECTED_SOURCES.append(tmpsource)
FAKE_SOURCES.append(tmpfakesource) FAKE_SOURCES.append(tmpfakesource)
immask.append(_immask)
matchmask.append(_matchmask)
bar.update() bar.update()
# To merge all the fake_sources and sources catalogs # To merge all the fake_sources and sources catalogs
...@@ -178,6 +196,9 @@ if nsources != 0: ...@@ -178,6 +196,9 @@ if nsources != 0:
_fake['find_peak'] = _fake['find_peak'] + n_detected _fake['find_peak'] = _fake['find_peak'] + n_detected
fake_sources = vstack([fake_sources, _fake]) fake_sources = vstack([fake_sources, _fake])
immask = np.sum(immask, axis=0)
print(np.unique(immask))
matchmask = np.sum(matchmask, axis=0)
fname = ('flux{:.2f}mJy_thresh{}_nsim{}.fits' fname = ('flux{:.2f}mJy_thresh{}_nsim{}.fits'
.format(_flux.to_value(unit=u.mJy), .format(_flux.to_value(unit=u.mJy),
...@@ -196,6 +217,8 @@ if nsources != 0: ...@@ -196,6 +217,8 @@ if nsources != 0:
hdul.append(fits.BinTableHDU(data=sources, hdul.append(fits.BinTableHDU(data=sources,
name='Detected_Sources')) name='Detected_Sources'))
hdul.append(fits.BinTableHDU(data=fake_sources, name='Fake_Sources')) hdul.append(fits.BinTableHDU(data=fake_sources, name='Fake_Sources'))
hdul.append(ImageHDU(data=immask, name='imagemask'))
hdul.append(ImageHDU(data=matchmask, name='matchmask'))
hdul = fits.HDUList(hdul) hdul = fits.HDUList(hdul)
hdul.writeto(outfile, overwrite=False) hdul.writeto(outfile, overwrite=False)
print('results written to {}'.format(outfile)) print('results written to {}'.format(outfile))
...@@ -215,14 +238,27 @@ if nsources == 0: ...@@ -215,14 +238,27 @@ if nsources == 0:
res = list(zip(*res)) res = list(zip(*res))
DETECTED_SOURCES = res[0] DETECTED_SOURCES = res[0]
immask = res[2]
matchmask = res[3]
else: else:
DETECTED_SOURCES = [] DETECTED_SOURCES = []
immask = []
matchmask = []
with ProgressBar(nsim) as bar: with ProgressBar(nsim) as bar:
for iloop in range(nsim): for iloop in range(nsim):
tmpsource, _ = helpfunc(jk_iter) (tmpsource,
tmpfakesource,
_immask,
_matchmask) = helpfunc(jk_iter)
DETECTED_SOURCES.append(tmpsource) DETECTED_SOURCES.append(tmpsource)
immask.append(_immask)
matchmask.append(_matchmask)
bar.update() bar.update()
immask = np.sum(immask, axis=0)
print(np.unique(immask))
matchmask = np.sum(matchmask, axis=0)
# To merge all the fake_sources and sources catalogs # To merge all the fake_sources and sources catalogs
sources = Table() sources = Table()
for _detected in DETECTED_SOURCES: for _detected in DETECTED_SOURCES:
...@@ -247,6 +283,10 @@ if nsources == 0: ...@@ -247,6 +283,10 @@ if nsources == 0:
if len(sources) > 0: if len(sources) > 0:
hdul.append(fits.BinTableHDU(data=sources, hdul.append(fits.BinTableHDU(data=sources,
name='Detected_Sources')) name='Detected_Sources'))
hdul.append(ImageHDU(data=immask, name='imagemask'))
hdul.append(ImageHDU(data=matchmask, name='matchmask'))
hdul = fits.HDUList(hdul) hdul = fits.HDUList(hdul)
hdul.writeto(outfile, overwrite=False) hdul.writeto(outfile, overwrite=False)
print('results written to {}'.format(outfile)) print('results written to {}'.format(outfile))
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