Commit 6e2d3f24 authored by LUSTIG Peter's avatar LUSTIG Peter

multiprocessing seed bug in catalog creation fixed

parent 71e22c18
......@@ -25,7 +25,10 @@ import sys
import datetime
from astropy.io.fits import ImageHDU
from copy import deepcopy, copy
# import dill
from utils import CombineMeasurements
#import dill as pickle
plt.close('all')
'''
%load_ext autoreload
......@@ -52,13 +55,14 @@ def get_size(obj, seen=None):
size += sum([get_size(k, seen) for k in obj.keys()])
elif hasattr(obj, '__dict__'):
size += get_size(obj.__dict__, seen)
elif hasattr(obj, '__iter__') and not isinstance(obj, (str, bytes, bytearray)):
elif hasattr(obj, '__iter__') and not isinstance(obj, (str, bytes,
bytearray)):
size += sum([get_size(i, seen) for i in obj])
return size
def fake_worker(img, 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,
retmask=True, photometry=False, **kwargs):
"""The completness purity worker, create a fake dataset from an jackknifed
image and return catalogs
......@@ -73,11 +77,8 @@ def fake_worker(img, min_threshold=2, nsources=8**2, flux=1*u.Jy,
arguments for source injection
"""
# print('using Jackknife object', jkiter)
'''
_jkiter = copy(jkiter)
print('using Jackknife object', _jkiter)
img = _jkiter(parity_threshold)
'''
np.random.seed(int(os.urandom(4).hex(), 16))
# img = IMAGE
# Renormalize the stddev
std = img.check_SNR()
......@@ -100,7 +101,7 @@ def fake_worker(img, min_threshold=2, nsources=8**2, flux=1*u.Jy,
# The gaussian fit from subpixel=True is very slow here...
mf_img.detect_sources(threshold=min_threshold)
if photometry:
if photometry and mf_img.sources is not None:
img.phot_sources(sources=mf_img.sources, peak=False, psf=True)
sources = img.sources
else:
......@@ -113,21 +114,16 @@ def fake_worker(img, min_threshold=2, nsources=8**2, flux=1*u.Jy,
else:
return sources, fake_sources
def CreateParallelJackkifes(jkiter, parity_threshold):
np.random.seed(int(os.urandom(4).hex(), 16))
jk = jkiter(parity_threshold=parity_threshold)
return jk
plt.close('all')
DATA_DIR_SERVER = Path("/data/NIKA/Reduced/"
"HLS091828_common_mode_one_block/v_1")
DATA_DIR_MYPC = Path("/home/peter/Dokumente/Uni/Paris/Stage/data/v_1")
flux = np.geomspace(1, 10, 3) * u.mJy
flux = [0]
flux = np.linspace(1, 3, 2)*u.mJy
##############################################################
flux = np.array([10]) * u.mJy
nsim = 4
nsim = 8
min_detection_threshold = 3
nsources = 2
#nsources = 0
photometry = True
ncores = 2
......@@ -136,15 +132,41 @@ parity_threshold = 1
outdir = Path("montecarlo_results/")
outdir = Path("testdir")
filenamecomment = ''
##############################################################
DATA_DIR_SERVER = Path("/data/NIKA/Reduced/"
"HLS091828_common_mode_one_block/v_1")
DATA_DIR_MYPC = Path("/home/peter/Dokumente/Uni/Paris/Stage/data/v_1")
message = '''\
###############################################################
# Running Monte Carlo Tests #
# #
# No Simulations per flux: {:5d} #
# Number of CPUs used: {:5d} #
# Minimum Detection Threshold:{:5.1f} #
# No Different Fluxes: {:5d} #
# Number of Injected Sources: {:5d} #
# Photometry: {} #
###############################################################
'''.format(nsim, ncores, min_detection_threshold, len(flux), nsources,
photometry)
print(message, "\n")
timeprefix = '{date:%y%m%d_%H%M%S}_'.format(date=datetime.datetime.now())
timeprefix += filenamecomment
if not outdir.exists():
print('creating directory {}'.format(outdir))
outdir.mkdir()
print("Creating Jackkife Object")
''' load maps for jackkifing'''
t0 = clock()
if DATA_DIR_SERVER.exists():
DATA_DIR = DATA_DIR_SERVER
elif DATA_DIR_MYPC.exists():
......@@ -152,142 +174,54 @@ elif DATA_DIR_MYPC.exists():
else:
sys.exit("Raw data path not found. Exit.")
jk_filenames = list(Path(DATA_DIR).glob('*/map.fits'))
for i in range(len(jk_filenames)):
jk_filenames[i] = str(jk_filenames[i])
message = '''\
###############################################################
# Running Monte Carlo Tests #
# #
# No Simulations per flux: {:5d} #
# Number of CPUs used: {:5d} #
# Minimum Detection Threshold:{:5.1f} #
# No Different Fluxes: {:5d} #
# Number of Injected Sources: {:5d} #
###############################################################
'''.format(nsim, ncores, min_detection_threshold, len(flux), nsources)
print(message, "\n")
print("Creating Jackkife Object")
t0 = clock()
jk_iter = Jackknife(jk_filenames, n=nsim)
print('Done in {:.2f}s'.format(clock()-t0))
print("Creating {} Jackkife Maps".format(nsim))
to = clock()
jackknifes = []
for i in range(nsim):
jackknifes.append(jk_iter(parity_threshold=parity_threshold))
print('Done in {:.2f}s'.format(clock()-t0))
min_source_dist = 2 * jackknifes[0].beam.fwhm_pix.value
p = Pool(ncores)
# set forced minimal distance of two fake sources
min_source_dist = 2 * jackknifes[0].beam.fwhm_pix.value
print('Begin Monte Carlo')
if nsources != 0:
for _flux in flux:
helpfunc = partial(fake_worker, **{'min_threshold':
min_detection_threshold,
'nsources': nsources, 'flux': _flux,
'within': (0, 1),
'cat_gen': pos_uniform,
'dist_threshold': min_source_dist,
'parity_threshold': 0.5,
'photometry': photometry})
print('Simulation with {:.2f}...'.format(_flux))
if(1):
res = p.map(helpfunc, jackknifes)
res = list(zip(*res))
(DETECTED_SOURCES,
FAKE_SOURCES,
immask,
matchmask) = res
else:
DETECTED_SOURCES = []
FAKE_SOURCES = []
immask = []
matchmask = []
with ProgressBar(nsim) as bar:
for iloop in range(nsim):
(tmpsource,
tmpfakesource,
_immask,
_matchmask) = helpfunc(jk_iter)
DETECTED_SOURCES.append(tmpsource)
FAKE_SOURCES.append(tmpfakesource)
immask.append(_immask)
matchmask.append(_matchmask)
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])
immask = np.sum(immask, axis=0)
print(np.unique(immask))
matchmask = np.sum(matchmask, axis=0)
fname = ('flux{:.2f}mJy_thresh{}_nsim{}.fits'
.format(_flux.to_value(unit=u.mJy),
min_detection_threshold, nsim))
fname = timeprefix + fname
outfile = outdir / fname
phdu = fits.PrimaryHDU(header=jackknifes[0].wcs.to_header())
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.append(ImageHDU(data=immask, name='imagemask'))
hdul.append(ImageHDU(data=matchmask, name='matchmask'))
hdul = fits.HDUList(hdul)
hdul.writeto(outfile, overwrite=False)
print('results written to {}'.format(outfile))
p = Pool(ncores)
if nsources == 0:
flux = np.array([0]) * u.mJy
for _flux in flux:
helpfunc = partial(fake_worker, **{'min_threshold':
min_detection_threshold,
'flux': None,
'nsources': 0,
'nsources': nsources, 'flux': _flux,
'within': (0, 1),
'parity_threshold': 0.5,
'cat_gen': pos_uniform,
'dist_threshold': min_source_dist,
'photometry': photometry})
print('Simulation without sources...')
if nsources != 0:
print('Simulation with {:.2f}...'.format(_flux))
else:
print('simulation without source injection')
if(1):
res = p.map(helpfunc, jackknifes)
res = list(zip(*res))
DETECTED_SOURCES = res[0]
immask = res[2]
matchmask = res[3]
(DETECTED_SOURCES,
FAKE_SOURCES,
immask,
matchmask) = res
else:
DETECTED_SOURCES = []
FAKE_SOURCES = []
immask = []
matchmask = []
with ProgressBar(nsim) as bar:
......@@ -295,43 +229,41 @@ if nsources == 0:
(tmpsource,
tmpfakesource,
_immask,
_matchmask) = helpfunc(jk_iter)
_matchmask) = helpfunc(jackknifes[iloop])
DETECTED_SOURCES.append(tmpsource)
FAKE_SOURCES.append(tmpfakesource)
immask.append(_immask)
matchmask.append(_matchmask)
bar.update()
# merge calalogues to one:
sources, fake_sources = CombineMeasurements(DETECTED_SOURCES, FAKE_SOURCES)
# create global masks:
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
sources = Table()
for _detected in DETECTED_SOURCES:
n_detected = len(sources)
if _detected is not None:
_detected['ID'] = _detected['ID'] + n_detected
sources = vstack([sources, _detected])
fname = ('nosources_thresh{}_nsim{}.fits'
.format(min_detection_threshold, nsim))
fname = ('flux{:.2f}mJy_thresh{}_nsources{}_nsim{}_phot{}.fits'
.format(_flux.to_value(unit=u.mJy),
min_detection_threshold, nsources, nsim, photometry))
fname = timeprefix + fname
outfile = outdir / fname
phdu = fits.PrimaryHDU(header=jackknifes[0].wcs.to_header())
phdu.header['influx'] = 0
phdu.header['influx'] = '{}'.format(_flux)
phdu.header['nsim'] = nsim
phdu.header['sourcespersim'] = 0
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'))
if len(fake_sources) > 0:
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.writeto(outfile, overwrite=False)
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