Commit e1cd14e9 authored by LUSTIG Peter's avatar LUSTIG Peter

Separated functions and script, added argument parser.

parent a6ea0ab6
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 18 11:08:28 2018
@author: peter
"""
from mpi4py import MPI
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from nikamap import NikaMap
#from completness import *
import myfunctions as mfct
from time import clock
import argparse
#import sys
import completness as cp
import os
plt.close('all')
flux_bins = 2
flux_range = [0.1, 10]
threshold_bins = 5
threshold_range = [2, 5]
threshold_range = [3, 5]
comm = MPI.COMM_WORLD
rank, size = comm.Get_rank(), comm.Get_size()
plt.ion()
desc = "Place false sources with a defined flux on Jackknife maps, detect sources and calculate the purity and completness of the detection."
parser = argparse.ArgumentParser(description=desc)
#parser.add_argument("-n", "--noisefile", action="store",
# dest="noisefile",
# help="Name of the Noise Containing Input File")
parser.add_argument(action="store",
dest="noisefile",
help="Name of the Noise Containing Input File")
parser.add_argument("-s", "--simulations", action="store",
dest="nsim", type=int, default=10,
help="Number of Simulations")
parser.add_argument("-f", "--flux", action="store",
dest="flux", type=float, default=10.,
help="Flux of the Created Sources in mJy")
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: completness-method-flux_[flux]-nsim_[nsim].fits")
parser.add_argument("-d", "--outdir", action="store",
dest="outdir", default="",
help="Output directory. Default: This folder.")
args=parser.parse_args()
nsim = args.nsim
flux = args.flux * u.mJy
noisefile = args.noisefile
method=args.method
outfile = args.outfile
outdir = args.outdir
#flux bins
#flux range
#threshold bins
#threshold range
############### Begin Calculation #############
sim_per_core = mfct.Distribute_Njobs(size, nsim)
data = NikaMap.read(noisefile)
fluxes = np.logspace(np.log10(flux_range[0]), np.log10(flux_range[1]), flux_bins)*u.mJy
fluxes_edges = np.logspace(np.log10(flux_range[0]), np.log10(flux_range[1]), flux_bins + 1)*u.mJy
threshold = np.linspace(threshold_range[0], threshold_range[1], threshold_bins)
threshold_edges = np.linspace(threshold_range[0], threshold_range[1], threshold_bins+1)
# Create 4D wcs
shape_4D, wcs_4D = cp.completness_purity_wcs(data.shape, data.wcs, bins=20,
flux_range=flux_range, flux_bins=flux_bins, flux_log=True,
threshold_range=threshold_range, threshold_bins=threshold_bins)
# Testing the lower edges
wcs_threshold = wcs_4D.sub([3])
assert np.all(np.abs(wcs_threshold.all_pix2world(np.arange(threshold_bins+1)-0.5, 0) - threshold_edges) < 1e-15)
wcs_flux = wcs_4D.sub([4])
assert np.all(np.abs(wcs_flux.all_pix2world(np.arange(flux_bins+1)-0.5, 0) - fluxes_edges.value) < 1e-13)
comm.Barrier()
if rank == 0:
T0 = clock()
comm.Barrier()
comp, comp_n, pur, pur_n = cp.completness_purity_2(flux, nsim=sim_per_core[rank],
nsources=8**2, within=(0, 1),
shape=shape_4D[0:3],
wcs=wcs_4D.sub([1, 2, 3]), jk_map=data)
comm.Barrier()
all_comp = comm.gather(comp, root=0)
all_comp_n = comm.gather(comp_n, root=0)
all_pur = comm.gather(pur, root=0)
all_pur_n = comm.gather(pur_n, root=0)
if rank == 0:
if (method == "unknown" and "/" in noisefile):
method = noisefile.split("/")[-2]
if outfile is None:
outfile = "completness-" + method + "-flux_{:04.0f}-nsim_{:03}.fits".format(flux.value, nsim)
if outdir != "":
if outdir[-1] != "/":
outdir = outdir + "/"
if not os.path.exists(outdir):
os.makedirs(outdir)
outfile = outdir + outfile
all_comp = np.sum(all_comp, axis=0)
all_comp_n = np.sum(all_comp_n, axis=0)
all_pur = np.sum(all_pur, axis=0)
all_pur_n = np.sum(all_pur_n, axis=0)
with np.errstate(divide='ignore', invalid='ignore'):
completness = all_comp / all_comp_n[..., np.newaxis]
purity = all_pur / all_pur_n
cp.WriteResult(outfile, threshold, completness, purity, nsim, wcs_4D, method, flux)
cp.Plot_CompPur(completness, purity, threshold, nsim=nsim, savename=outfile[0:-5]+"-plot.png")
T = clock()-T0
print('Total calculation time: {:.0f}'.format(T))
comm.Barrier()
......@@ -33,10 +33,7 @@ if '__IPYTHON__' in globals():
ipython.magic('matplotlib tk')
comm = MPI.COMM_WORLD
rank, size = comm.Get_rank(), comm.Get_size()
plt.ion()
DATA_DIR = "/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/ForMarseille/HLS091828_common_mode_one_block/"
......@@ -240,8 +237,8 @@ def completness_purity(flux, wcs=None, shape=None, nsources=8**2, within=(1 / 4,
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)
#FAKE_SOURCES.append(fake_sources)
#DETECTED_SOURCES.append(sources)
_completness, _norm_comp = completness_worker(shape, wcs, sources, fake_sources, min_threshold, max_threshold)
......@@ -286,8 +283,8 @@ def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, with
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)
# Save global array (MB requests)
FAKE_SOURCES.append(fake_sources)
DETECTED_SOURCES.append(sources)
#FAKE_SOURCES.append(fake_sources)
#DETECTED_SOURCES.append(sources)
_completness, _norm_comp = completness_worker(shape, wcs, sources, fake_sources, min_threshold, max_threshold)
......@@ -309,7 +306,7 @@ def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, with
def Plot_CompPur(completness, purity, threshold,nsim=None):
def Plot_CompPur(completness, purity, threshold,nsim=None, savename=None):
threshold_bins = completness.shape[-1]
fig, axes = plt.subplots(nrows=2, ncols=threshold_bins, sharex=True, sharey=True)
for i in range(threshold_bins):
......@@ -320,26 +317,28 @@ def Plot_CompPur(completness, purity, threshold,nsim=None):
axes[0, 0].set_title("{} simulations".format(nsim))
axes[0, 0].set_ylabel("completness")
axes[1, 0].set_ylabel("purity")
if savename is not None:
plt.savefig(savename)
plt.show(block=True)
def AddParametersToHeader(header, threshold, nsim, method=None):
def AddParametersToHeader(header, threshold, nsim, method, flux):
header['thresh'] = threshold
header['nsim'] = nsim
if method is not None:
header['method'] = method
header['method'] = method
header['flux_mJy'] = flux.value
def WriteResult(filename, threshold, completness, purity, nsim, wcs4D):
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)
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)
AddParametersToHeader(hdupur.header, thresh, nsim, method, flux)
allhdupur.append(hdupur)
wcs2D = wcs4D.sub([1,2])
......@@ -351,109 +350,6 @@ def WriteResult(filename, threshold, completness, purity, nsim, wcs4D):
hdul.writeto(filename, overwrite=True)
plt.close('all')
filenames = list(Path(DATA_DIR).glob('../*/map.fits'))
#filenames = DATA_DIR+'map.fits'
data = NikaMap.read(Path(DATA_DIR) / 'map_JK.fits')
#
# _data = next(Jackknife(filenames, n=None))
# # TODO: Should in principle be the same, but is not... check....
# _ = plt.hist((data.data - _data.data)[~data.mask & ~_data.mask], bins=1000, range=[-0.1, 0.1], log=True)
# Create the flux and threshold axes...
flux_bins = 2
flux_range = [0.1, 10]
# Does not really make sense... better define edges
fluxes = np.logspace(np.log10(flux_range[0]), np.log10(flux_range[1]), flux_bins)*u.mJy
fluxes_edges = np.logspace(np.log10(flux_range[0]), np.log10(flux_range[1]), flux_bins + 1)*u.mJy
# Does not really make sense... better define edges
threshold_bins = 5
threshold_range = [2, 5]
threshold_range = [3, 5]
threshold = np.linspace(threshold_range[0], threshold_range[1], threshold_bins)
threshold_edges = np.linspace(threshold_range[0], threshold_range[1], threshold_bins+1)
shape_4D, wcs_4D = completness_purity_wcs(data.shape, data.wcs, bins=20,
flux_range=flux_range, flux_bins=flux_bins, flux_log=True,
threshold_range=threshold_range, threshold_bins=threshold_bins)
#print(wcs_4D)
# Testing the lower edges
wcs_threshold = wcs_4D.sub([3])
assert np.all(np.abs(wcs_threshold.all_pix2world(np.arange(threshold_bins+1)-0.5, 0) - threshold_edges) < 1e-15)
wcs_flux = wcs_4D.sub([4])
assert np.all(np.abs(wcs_flux.all_pix2world(np.arange(flux_bins+1)-0.5, 0) - fluxes_edges.value) < 1e-13)
FAKE_SOURCES = []
DETECTED_SOURCES = []
#wcsheader = wcs_4D.to_header()
#hdu = fits.PrimaryHDU(header=wcsheader)
#hdu.writeto('wcs.fits')
#print('written')
# DEBUG :
# flux, nsources, within, wcs, shape, nsim, jk_filenames = 10*u.mJy, 8**2, (0, 1), wcs_4D.sub([1, 2, 3]), shape_4D[0:3], np.multiply(*shape_4D[0:2]) * 100, filenames
# This is a single run check for a single flux
nsim = 4
#ncores = 3
allflux = [10*u.mJy]
sim_per_core = mfct.Distribute_Njobs(size, nsim)
print(sim_per_core[rank])
#completness_purity_2 = partial(completness_purity_2, **{'nsources':8**2, 'within':(0, 1),
# 'wcs':wcs_4D.sub([1, 2, 3]),
# 'shape':shape_4D[0:3],
# 'jk_map':data,
# 'nsim':2})
comm.Barrier()
if rank == 0:
T0 = clock()
for irun, flux in enumerate(allflux):
comm.Barrier()
if rank ==0:
t0 = clock()
comp, comp_n, pur, pur_n = completness_purity_2(10*u.mJy, nsim=sim_per_core[rank],
nsources=8**2, within=(0, 1),
shape=shape_4D[0:3],
wcs=wcs_4D.sub([1, 2, 3]), jk_map=data)
comm.Barrier()
all_comp = comm.gather(comp, root=0)
all_comp_n = comm.gather(comp_n, root=0)
all_pur = comm.gather(pur, root=0)
all_pur_n = comm.gather(pur_n, root=0)
if rank == 0:
all_comp = np.sum(all_comp, axis=0)
all_comp_n = np.sum(all_comp_n, axis=0)
all_pur = np.sum(all_pur, axis=0)
all_pur_n = np.sum(all_pur_n, axis=0)
with np.errstate(divide='ignore', invalid='ignore'):
completness = all_comp / all_comp_n[..., np.newaxis]
purity = all_pur / all_pur_n
WriteResult('test.fits', threshold, completness, purity, nsim, wcs_4D)
Plot_CompPur(completness, purity, threshold, nsim=nsim)
# Save result !!
......
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