Commit 33497262 authored by ljbing's avatar ljbing
Browse files

for PIIC outputs

parent e16ed33a
from __future__ import absolute_import, division, print_function
from pathlib import Path
import os,pdb
import uuid
import numpy as np
import matplotlib.pyplot as plt
from multiprocess 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
from time import process_time
import argparse
import sys
import datetime
from astropy.io.fits import ImageHDU
from copy import deepcopy, copy
from utils import CombineMeasurements
import binascii
# import dill as pickle
plt.close("all")
"""
%load_ext autoreload
%autoreload 2
%matplotlib tk
"""
plt.ion()
def worker(
img,
min_threshold=2,
nsources=8 ** 2,
flux=1 * u.Jy,
within=(0, 1),
cat_gen=pos_uniform,
retmask=True,
detection=True,
photometry=False,
**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
"""
# print('using Jackknife object', jkiter)
assert detection or photometry, "you have to chose either detection or " "photometry or both"
# np.random.seed(int(os.urandom(4).hex(), 16))
np.random.seed(int(binascii.hexlify(os.urandom(4)), 16))
# img = IMAGE
# Renormalize the stddev
std = img.check_SNR()
img.uncertainty.array *= std
# Actually rather slow... maybe check the code ?
# print(flux)
if flux is not None and (nsources != 0):
img.add_gaussian_sources(nsources=nsources, within=within, peak_flux=flux, cat_gen=cat_gen, **kwargs)
if detection:
# ... 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)
if photometry:
if detection and mf_img.sources is not None:
img.phot_sources(sources=mf_img.sources, peak=False, psf=True)# changed to get both peak and psf fluxes
mf_img.phot_sources(sources=img.sources, peak=True, psf=False)
sources = mf_img.sources
if not detection:
# print('testprint', img.sources.keys())
img.phot_sources(sources=img.fake_sources, peak=False, psf=True)
sources = img.sources
else:
sources = mf_img.sources
if detection:
fake_sources = mf_img.fake_sources
else:
fake_sources = img.fake_sources
if detection:
mf_img_mask = mf_img.mask
else:
mf_img_mask = np.zeros(img.data.shape)
if retmask:
return sources, fake_sources, img.mask, mf_img_mask
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
##############################################################
def main(
data_dir,
output_dir=None,
pattern="**/*.fits",
label=None,
min_detection_threshold=4,
photometry=True,
detection=True,
band="1mm",
):
"""Process all maps in a given directory
Parameters
----------
data_dir : Path or str
the input path for the jackknifes maps
output_dir : Path or str, optional
the output directory, by default '.'
pattern : str, optional
the pattern to look for maps, by default '**/*.fits', all fits file in all sub directories
label : str
the output label for this jackknife
min_detection_threshold : float
minimum detection thresdhold, by default 4
photometry : bool
do the photometry of the sources (?), by default True
detection : bool
do the detection of the sources (?), by default True
band : str, optionnal
the band to real (for IDL maps), by default '1mm'
"""
if label is None:
label = str(uuid.uuid4())
if output_dir is None:
output_dir = Path(".")
output_dir = Path(output_dir) / f'{label}_3' / f'{band}'
if not output_dir.exists():
print("creating directory {}".format(output_dir))
os.makedirs(output_dir)
# TODO: Shall we pass that as argument also ?
fluxes = np.geomspace(0.05, 5, 50) * u.mJy # np.geomspace(0.1, 10, 50) for 1mm
nsim = 2000
nsources = 50
parity_threshold = 0 # incomplete map when using 1 as threshold
filenamecomment = "NewRealisation_"
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: {} #
# Detection: {} #
# Band: {} #
###############################################################
""".format(
nsim, cpu_count(), min_detection_threshold, len(fluxes), nsources, photometry, detection, band
)
print(message, "\n")
timeprefix = "{date:%y%m%d_%H%M%S}_".format(date=datetime.datetime.now())
timeprefix += filenamecomment
print("Creating Jackkife Object")
t0 = process_time()
jk_filenames = Path(data_dir).glob("*/map.fits")
jk_iter = Jackknife(jk_filenames, band=band, n=nsim, parity_threshold=parity_threshold)
print("Done in {:.2f}s".format(process_time() - t0))
print("Begin Monte Carlo")
if nsources == 0:
fluxes = np.array([0]) * u.mJy
print("Creating {} Jackkife Maps".format(nsim))
to = process_time()
jackknifes = []
jk_std = []
# TODO: This take a LOT of time and could be done in parallel in the worker, as needed, need to see how to pass the Jackknife object to the worker...
for i in ProgressBar(range(nsim)):
jackknifes.append(jk_iter())
jk_std.append(jackknifes[i].data)
jk_std = np.std(jk_std,axis=0)
print("Done in {:.2f}s".format(process_time() - t0))
#pdb.set_trace()
for flux in fluxes:
#if flux.value < 0.098:
# print(f'skip {flux}')
# continue
min_source_dist = 2 * jackknifes[0].beam.fwhm_pix.value
helpfunc = partial(
worker,
**{
"min_threshold": min_detection_threshold,
"nsources": nsources,
"flux": flux,
"within": (0, 1),
"cat_gen": pos_uniform,
"dist_threshold": min_source_dist,
"photometry": photometry,
"detection": detection,
},
)
if nsources != 0:
print("Simulation with {:.2f}...".format(flux))
else:
print("simulation without source injection")
with Pool(cpu_count()) as p:
res = p.map(helpfunc, jackknifes) # func, args
#pdb.set_trace()
# Should work, but don't... investigate :
# res = ProgressBar.map(helpfunc, jackknifes, multiprocess=cpu_count())
res = list(zip(*res))
(DETECTED_SOURCES, FAKE_SOURCES, immask, matchmask) = res
# merge calalogues to one:
sources, fake_sources = CombineMeasurements(DETECTED_SOURCES, FAKE_SOURCES)
# create global masks:
immask = np.sum(immask, axis=0)
matchmask = np.sum(matchmask, axis=0)
fname = "flux{:.3f}mJy_thresh{}_nsources{}_nsim{}_phot{}_detect{}.fits".format(
flux.to_value(unit=u.mJy), min_detection_threshold, nsources, nsim, photometry, detection
)
fname = timeprefix + fname
outfile = output_dir / fname
stdname = f"stddev.fits"
stdfile = output_dir / stdname
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"))
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))
stdhdu = fits.PrimaryHDU()
stdhdu = [stdhdu]
stdhdu.append(ImageHDU(data=jk_std, name="stddev_from_jackknife", header=jackknifes[0].wcs.to_header()))
stdhdu = fits.HDUList(stdhdu)
stdhdu.writeto(stdfile, overwrite=True)
if __name__ == "__main__":
##############################################################
# DATA_DIR_SERVER = Path("/data/NIKA/Reduced/"
# "HLS091828_common_mode_one_block/v_1")
# DATA_DIR_SERVER = Path("/data/NIKA/Reduced/May2018/"
# "HLS091828_common_mode_one_block/v_1/")
# DATA_DIR_SERVER = Path("/data/NIKA/Reduced/May2018/HLS091828_atm_and"
# "_all_box_iter_IterativeMM_CosSinElevOffset/v_1")
# DATA_DIR_MYPC = Path("/home/peter/Dokumente/Uni/Paris/Stage/data/v_1")
band = "2mm"
method = 211
method_dir = Path(f"{method}")
min_sn = 3.
DATA_DIR_SERVER = Path(f"/data/PHOTOM/Reduced/23515M/GOODSNORTH/{method}/iter15/v_1")
DATA_DIR_MYPC = Path("/home/lbing/data/purity_completeness")
outdir = DATA_DIR_MYPC
main(DATA_DIR_SERVER, output_dir=outdir, label=str(method), band=band,min_detection_threshold=min_sn)
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