Commit 221a157d authored by LUSTIG Peter's avatar LUSTIG Peter

Catalog of detected sources is written in output file

parent ce03c1cc
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import dill as pickle
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
......@@ -11,6 +12,7 @@ from time import clock
from multiprocess import Pool
from functools import partial
import sys
from astropy.table import vstack
#import sys
import completness as cp
......@@ -25,14 +27,14 @@ threshold_bins = 5
threshold_range = [2, 10]
#threshold_range = [3, 5]
nsim = 6
nsim = 10
flux = 10 * u.mJy
noisefile = "../../N2CLS/ForMarseille/HLS091828_common_mode_kids_out/map_JK.fits"
method="unknown"
outfile = None
outdir = "results/"
ncore = 3
ncore = 4
'''
import argparse
......@@ -118,6 +120,7 @@ all_comp = np.sum(tmpres[0], axis=0)
all_comp_n = np.sum(tmpres[1], axis=0)
all_pur = np.sum(tmpres[2], axis=0)
all_pur_n = np.sum(tmpres[3], axis=0)
sourceresult = vstack(tmpres[4])
print(all_comp.shape)
print(all_comp_n.shape)
......@@ -143,7 +146,7 @@ if outdir != "":
outfile = outdir + outfile
cp.WriteResult(outfile, threshold, completness, purity, nsim, wcs_4D, method, flux)
mfct.WriteResult(outfile, threshold, completness, purity, nsim, wcs_4D, method, flux, sourceresult)
cp.Plot_CompPur(completness, purity, threshold, nsim=nsim, savename=outfile[0:-5]+"-plot.png")
T = clock()-T0
......
......@@ -5,14 +5,14 @@ import os
import numpy as np
import matplotlib.pyplot as plt
from multiprocess import Pool, cpu_count
#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
from astropy.table import Table, vstack
from scipy.optimize import curve_fit
......@@ -84,7 +84,7 @@ def completness_purity_wcs(shape, wcs, bins=30,
return (bins, bins, threshold_bins, flux_bins), WCS(header)
def fake_worker(img, min_threshold=2, nsources=8**2, flux=1*u.Jy, within=(1 / 4, 3 / 4), pos_gen=pos_uniform, **kwargs):
def fake_worker(img, min_threshold=2, nsources=8**2, flux=1*u.Jy, within=(1 / 4, 3 / 4), pos_gen=pos_uniform, photmethod='peak',**kwargs):
"""The completness purity worker, create a fake dataset from an jackknifed image, detect sources
in doped image and return catalogs
......@@ -118,9 +118,27 @@ def fake_worker(img, min_threshold=2, nsources=8**2, flux=1*u.Jy, within=(1 / 4,
# ... and detect sources with the lowest threshold...
# The gaussian fit from subpixel=True is very slow here...
mf_img.detect_sources(threshold=min_threshold)
FluxResolution(mf_img, method=photmethod)
mf_img.phot_sources(peak=True, psf=False)
#flux = FluxResolution(mf_img)
return mf_img.sources, mf_img.fake_sources
def FluxResolution(nm, method='peak'):
if method== 'peak':
peak, psf = True, False
else:
peak, psf = False, True
nm.phot_sources(sources=nm.sources, peak=peak, psf=psf)
#np.savetxt('fluxhisttest.txt', nm.sources['flux_'+method])
return 0
def completness_worker(shape, wcs, sources, fake_sources, min_threshold=2, max_threshold=5):
"""Compute completness from the fake source catalog
......@@ -208,52 +226,53 @@ def purity_worker(shape, wcs, sources, max_threshold=2):
return _purity, _norm_pur
def completness_purity(flux, wcs=None, shape=None, nsources=8**2, within=(1 / 4, 3 / 4), nsim=500, jk_filenames=None):
"""Compute completness map for a given flux"""
print(flux)
# wcs_celestial = wcs.celestial
# Lower and upper edges ... take the center of the pixel for the upper edge
min_threshold, max_threshold = wcs.sub([3]).all_pix2world([-0.5, shape[2]-1], 0)[0]
completness = np.zeros(shape, dtype=np.float)
norm_comp = np.zeros(shape[0:2], dtype=np.float)
purity = np.zeros(shape, dtype=np.float)
norm_pur = np.zeros(shape, dtype=np.float)
jk_iter = Jackknife(jk_filenames, n=nsim)
print("Begin loop")
for img in ProgressBar(jk_iter):
print("In loop")
# %load_ext snakeviz
# %snakeviz the following line.... all is spend in the find_peaks / fit_2d_gaussian
# TODO: Change the find_peaks routine, or maybe just the fit_2d_gaussian to be FAST ! (Maybe look into gcntrd.pro routine or photutils.centroid.centroid_1dg maybe ?)
sources, fake_sources = fake_worker(img, nsources=nsources, within=within, flux=flux, pos_gen=pos_uniform, min_threshold=min_threshold, dist_threshold=img.beam.fwhm_pix.value)
# Save global array (MB requests)
#FAKE_SOURCES.append(fake_sources)
#DETECTED_SOURCES.append(sources)
_completness, _norm_comp = completness_worker(shape, wcs, sources, fake_sources, min_threshold, max_threshold)
completness += _completness
norm_comp += _norm_comp
_purity, _norm_pur = purity_worker(shape, wcs, sources, max_threshold)
purity += _purity
norm_pur += _norm_pur
# norm can be 0, so to avoid warning on invalid values...
with np.errstate(divide='ignore', invalid='ignore'):
completness /= norm_comp[..., np.newaxis]
purity /= norm_pur
# TODO: One should probably return completness AND norm if one want to combine several fluxes
return completness, purity
#def completness_purity(flux, wcs=None, shape=None, nsources=8**2, within=(1 / 4, 3 / 4), nsim=500, jk_filenames=None):
# """Compute completness map for a given flux"""
#
# print(flux)
# #fluxlist=[]
#
# # wcs_celestial = wcs.celestial
# # Lower and upper edges ... take the center of the pixel for the upper edge
# min_threshold, max_threshold = wcs.sub([3]).all_pix2world([-0.5, shape[2]-1], 0)[0]
#
# completness = np.zeros(shape, dtype=np.float)
# norm_comp = np.zeros(shape[0:2], dtype=np.float)
#
# purity = np.zeros(shape, dtype=np.float)
# norm_pur = np.zeros(shape, dtype=np.float)
#
# jk_iter = Jackknife(jk_filenames, n=nsim)
# print("Begin loop")
#
# for img in ProgressBar(jk_iter):
# print("In loop")
# # %load_ext snakeviz
# # %snakeviz the following line.... all is spend in the find_peaks / fit_2d_gaussian
# # TODO: Change the find_peaks routine, or maybe just the fit_2d_gaussian to be FAST ! (Maybe look into gcntrd.pro routine or photutils.centroid.centroid_1dg maybe ?)
# sources, fake_sources = fake_worker(img, nsources=nsources, within=within, flux=flux, pos_gen=pos_uniform, min_threshold=min_threshold, dist_threshold=img.beam.fwhm_pix.value)
#
# # Save global array (MB requests)
# #FAKE_SOURCES.append(fake_sources)
# #DETECTED_SOURCES.append(sources)
#
# _completness, _norm_comp = completness_worker(shape, wcs, sources, fake_sources, min_threshold, max_threshold)
#
# completness += _completness
# norm_comp += _norm_comp
#
# _purity, _norm_pur = purity_worker(shape, wcs, sources, max_threshold)
#
# purity += _purity
# norm_pur += _norm_pur
#
# # norm can be 0, so to avoid warning on invalid values...
# with np.errstate(divide='ignore', invalid='ignore'):
# completness /= norm_comp[..., np.newaxis]
# purity /= norm_pur
#
# # TODO: One should probably return completness AND norm if one want to combine several fluxes
# return completness, purity
def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, within=(1 / 4, 3 / 4), jk_map=None):
......@@ -271,6 +290,7 @@ def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, with
purity = np.zeros(shape, dtype=np.float)
norm_pur = np.zeros(shape, dtype=np.float)
sourcetable = Table()
#jk_iter = Jackknife(jk_filenames, n=nsim)
......@@ -279,6 +299,7 @@ def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, with
# %snakeviz the following line.... all is spend in the find_peaks / fit_2d_gaussian
# TODO: Change the find_peaks routine, or maybe just the fit_2d_gaussian to be FAST ! (Maybe look into gcntrd.pro routine or photutils.centroid.centroid_1dg maybe ?)
sources, fake_sources = fake_worker(jk_map, nsources=nsources, within=within, flux=flux, pos_gen=pos_uniform, min_threshold=min_threshold, dist_threshold=jk_map.beam.fwhm_pix.value)
sourcetable = vstack([sourcetable,sources])
# Save global array (MB requests)
#FAKE_SOURCES.append(fake_sources)
......@@ -293,6 +314,8 @@ def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, with
purity += _purity
norm_pur += _norm_pur
# norm can be 0, so to avoid warning on invalid values...
#with np.errstate(divide='ignore', invalid='ignore'):
......@@ -300,7 +323,7 @@ def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, with
# purity /= norm_pur
# TODO: One should probably return completness AND norm if one want to combine several fluxes
return completness, norm_comp, purity, norm_pur
return completness, norm_comp, purity, norm_pur, sourcetable
......@@ -326,33 +349,6 @@ def Plot_CompPur(completness, purity, threshold,nsim=None, savename=None):
#plt.show(block=True)
def AddParametersToHeader(header, threshold, nsim, method, flux):
header['thresh'] = threshold
header['nsim'] = nsim
header['method'] = method
header['flux_mJy'] = flux.value
def WriteResult(filename, threshold, completness, purity, nsim, wcs4D, method, flux):
allhducomp, allhdupur = [],[]
for i, thresh in enumerate(threshold):
hducomp = fits.ImageHDU(data=completness[i], name='completness_thresh_{:.1f}'.format(float(threshold[i])))
AddParametersToHeader(hducomp.header, thresh, nsim, method, flux)
allhducomp.append(hducomp)
hdupur=fits.ImageHDU(data=purity[i], name='purity_thresh_{:.1f}'.format(float(threshold[i])))
AddParametersToHeader(hdupur.header, thresh, nsim, method, flux)
allhdupur.append(hdupur)
wcs2D = wcs4D.sub([1,2])
hdu = fits.PrimaryHDU(header=wcs2D.to_header())
hdul = [hdu]+allhducomp+allhdupur
hdul = fits.HDUList(hdul)
hdul.writeto(filename, overwrite=True)
# Save result !!
......
......@@ -2,9 +2,37 @@
# -*- coding: utf-8 -*-
import numpy as np
from astropy.io import fits
def Distribute_Njobs(n_cpu, n_jobs):
n = np.repeat(int(n_jobs/n_cpu), n_cpu)
n[0:n_jobs%n_cpu] += 1
return n
def AddParametersToHeader(header, threshold, nsim, method, flux):
header['thresh'] = threshold
header['nsim'] = nsim
header['method'] = method
header['flux_mJy'] = flux.value
def WriteResult(filename, threshold, completness, purity, nsim, wcs4D, method, flux, sourcetable):
allhducomp, allhdupur = [],[]
for i, thresh in enumerate(threshold):
hducomp = fits.ImageHDU(data=completness[i], name='completness_thresh_{:.1f}'.format(float(threshold[i])))
AddParametersToHeader(hducomp.header, thresh, nsim, method, flux)
allhducomp.append(hducomp)
hdupur=fits.ImageHDU(data=purity[i], name='purity_thresh_{:.1f}'.format(float(threshold[i])))
AddParametersToHeader(hdupur.header, thresh, nsim, method, flux)
allhdupur.append(hdupur)
wcs2D = wcs4D.sub([1,2])
hdu = fits.PrimaryHDU(header=wcs2D.to_header())
hdul = [hdu]+allhducomp+allhdupur+[fits.BinTableHDU(sourcetable)]
hdul = fits.HDUList(hdul)
hdul.writeto(filename, 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