Commit 6d221375 authored by tahumada's avatar tahumada

refined levels and frequency changes

parent 1cd1a9b0
Pipeline #95883 failed with stages
in 39 minutes and 52 seconds
......@@ -292,7 +292,7 @@ def localize(
min_inclination=0, max_inclination=np.pi / 2,
min_distance=None, max_distance=None, prior_distance_power=None,
cosmology=False, mcmc=False, chain_dump=None,
enable_snr_series=True, f_high_truncate=0.95):
enable_snr_series=True, f_high_truncate=0.95,refine_level=1):
"""Localize a compact binary signal using the BAYESTAR algorithm.
Parameters
......@@ -321,7 +321,9 @@ def localize(
Truncate the noise power spectral densities at this factor times the
highest sampled frequency to suppress artifacts caused by incorrect
PSD conditioning by some matched filter pipelines.
refine_level: int, optional
1 or 0. If 1 will loop over 11 levels to create a Multi resolution fits file.
(default: 1, will preform the loop).
Returns
-------
skymap : `astropy.table.Table`
......@@ -368,7 +370,7 @@ def localize(
else:
args = (min_inclination, max_inclination, min_distance, max_distance,
prior_distance_power, cosmology, gmst, sample_rate, toas,
snr_series, responses, locations, horizons)
snr_series, responses, locations, horizons, refine_level)
skymap, log_bci, log_bsn = core.toa_phoa_snr(*args)
skymap = Table(skymap)
skymap.meta['log_bci'] = log_bci
......
......@@ -244,7 +244,10 @@ def sngl_inspiral_psd(waveform, mass1, mass2, f_min=10, f_max=2048, f_ref=0,
waveform = 'SEOBNRv4_ROM'
approx, ampo, phaseo = get_approximant_and_orders_from_string(waveform)
log.info('Selected template: %s', waveform)
## control
print('#######################\n FREQUENCY',float(kwargs.get('f_final') or f_max))
print(kwargs)
# Generate conditioned template.
params = lal.CreateDict()
lalsimulation.SimInspiralWaveformParamsInsertPNPhaseOrder(params, phaseo)
......@@ -259,7 +262,9 @@ def sngl_inspiral_psd(waveform, mass1, mass2, f_min=10, f_max=2048, f_ref=0,
S2z=float(kwargs.get('spin2z') or 0),
distance=1e6 * lal.PC_SI, inclination=0, phiRef=0,
longAscNodes=0, eccentricity=0, meanPerAno=0,
deltaF=0, f_min=f_min, f_max=f_max, f_ref=f_ref,
deltaF=0, f_min=f_min,
f_max=float(kwargs.get('f_final') or f_max),
f_ref=f_ref,
LALparams=params, approximant=approx)
# Force `plus' and `cross' waveform to be in quadrature.
......
......@@ -85,7 +85,7 @@ class LigoLWEventSource(OrderedDict, EventSource):
self._make_events(doc, psd_file, coinc_def))
_template_keys = '''mass1 mass2
spin1x spin1y spin1z spin2x spin2y spin2z'''.split()
spin1x spin1y spin1z spin2x spin2y spin2z f_final'''.split()
_invert_phases = {
'pycbc': False,
......@@ -193,6 +193,10 @@ class LigoLWEventSource(OrderedDict, EventSource):
raise ValueError(
'Template arguments are not identical for all detectors!')
template_args = template_args[0]
# control
print('######\n')#,template_args)
for i,arg in enumerate(template_args):
print(self._template_keys[i],template_args[arg])
invert_phases = self._phase_convention(
program_for_process_id[coinc.process_id])
......
......@@ -67,6 +67,9 @@ def parser():
parser.add_argument(
'--condor-submit', action='store_true',
help='submit to Condor instead of running locally')
parser.add_argument(
'--refine-level', type=int,required=False,
help='boolean (0 or 1) for produce a multi order fits file', default=1)
return parser
......@@ -100,7 +103,7 @@ def main(args=None):
event_source = events.open(*opts.input, sample=opts.pycbc_sample)
mkpath(opts.output)
if opts.condor_submit:
if opts.seed is not None:
raise NotImplementedError(
......@@ -161,7 +164,9 @@ def main(args=None):
opts.max_distance, opts.prior_distance_power,
opts.cosmology, mcmc=opts.mcmc, chain_dump=chain_dump,
enable_snr_series=opts.enable_snr_series,
f_high_truncate=opts.f_high_truncate)
f_high_truncate=opts.f_high_truncate,
refine_level=opts.refine_level
)
sky_map.meta['objid'] = coinc_event_id
except (ArithmeticError, ValueError):
log.exception('%s:sky localization failed', coinc_event_id)
......
......@@ -47,7 +47,7 @@ def parser():
'--df', metavar='Hz', type=float, default=1.0,
help='Frequency step size [default: %(default)s]')
parser.add_argument(
'--f-max', metavar='Hz', type=float, default=2048.0,
'--f-max', metavar='Hz', type=float, default=2048.0,#default=512.0,
help='Maximum frequency [default: %(default)s]')
detector_group = parser.add_argument_group(
......
......@@ -854,7 +854,8 @@ bayestar_pixel *bayestar_sky_map_toa_phoa_snr(
const float complex **snrs, /* Complex SNR series */
const float (**responses)[3], /* Detector responses */
const double **locations, /* Barycentered Cartesian geographic detector positions (light seconds) */
const double *horizons /* SNR=1 horizon distances for each detector */
const double *horizons, /* SNR=1 horizon distances for each detector */
int refine_level /* Bool to go into multi order loop, default = 1 */
) {
/* Initialize precalculated tables. */
bayestar_init();
......@@ -1058,36 +1059,42 @@ bayestar_pixel *bayestar_sky_map_toa_phoa_snr(
/* Sort pixels by ascending posterior probability. */
bayestar_pixels_sort_prob(pixels, len);
/* Adaptively refine until order=11 (nside=2048). */
for (unsigned char level = order0; level < 11; level ++)
/* int refine_level = 0; */
/*Checking for loop access*/
printf("########################## \n %d \n",refine_level);
if (refine_level)
{
/* Adaptively refine the pixels that contain the most probability. */
pixels = bayestar_pixels_refine(pixels, &len, npix0 / 4);
if (!pixels)
goto done;
ITT_TASK_BEGIN(itt_domain, itt_task_refinement_step);
#pragma omp parallel for schedule(guided)
for (unsigned long i = len - npix0; i < len; i ++)
/* Adaptively refine until order=11 (nside=2048). */
for (unsigned char level = order0; level < 6; level ++)
{
/* Adaptively refine the pixels that contain the most probability. */
printf("#### doing refinement at level %d\n", level);
pixels = bayestar_pixels_refine(pixels, &len, npix0 / 4);
if (!pixels)
goto done;
ITT_TASK_BEGIN(itt_domain, itt_task_refinement_step);
#pragma omp parallel for schedule(guided)
for (unsigned long i = len - npix0; i < len; i ++)
{
if (OMP_WAS_INTERRUPTED)
OMP_EXIT_LOOP_EARLY;
bayestar_sky_map_toa_phoa_snr_pixel(integrators, 1, pixels[i].uniq,
pixels[i].value, gmst, nifos, nsamples, n_u_points,
sample_rate, epochs, snrs, responses, locations, horizons,
u_points_weights);
}
ITT_TASK_END(itt_domain);
if (OMP_WAS_INTERRUPTED)
OMP_EXIT_LOOP_EARLY;
goto done;
bayestar_sky_map_toa_phoa_snr_pixel(integrators, 1, pixels[i].uniq,
pixels[i].value, gmst, nifos, nsamples, n_u_points,
sample_rate, epochs, snrs, responses, locations, horizons,
u_points_weights);
/* Sort pixels by ascending posterior probability. */
bayestar_pixels_sort_prob(pixels, len);
}
ITT_TASK_END(itt_domain);
if (OMP_WAS_INTERRUPTED)
goto done;
/* Sort pixels by ascending posterior probability. */
bayestar_pixels_sort_prob(pixels, len);
}
/* Evaluate distance layers. */
ITT_TASK_BEGIN(itt_domain, itt_task_final_step);
#pragma omp parallel for schedule(guided)
......
......@@ -102,7 +102,8 @@ bayestar_pixel *bayestar_sky_map_toa_phoa_snr(
const float complex **snrs, /* Complex SNR series */
const float (**responses)[3], /* Detector responses */
const double **locations, /* Barycentered Cartesian geographic detector positions (light seconds) */
const double *horizons /* SNR=1 horizon distances for each detector */
const double *horizons, /* SNR=1 horizon distances for each detector */
int refine_level /* Bool to go into multi order loop, default = 1 */
);
double bayestar_log_posterior_toa_phoa_snr(
......
......@@ -721,21 +721,22 @@ static PyObject *sky_map_toa_phoa_snr(
PyObject *responses_obj;
PyObject *locations_obj;
PyObject *horizons_obj;
int refine_level;
/* Names of arguments */
static const char *keywords[] = {"min_inclination", "max_inclination",
"min_distance", "max_distance", "prior_distance_power", "cosmology",
"gmst", "sample_rate", "epochs", "snrs", "responses", "locations",
"horizons", NULL};
"horizons", "refine_level", NULL};
/* Parse arguments */
/* FIXME: PyArg_ParseTupleAndKeywords should expect keywords to be const */
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wincompatible-pointer-types"
if (!PyArg_ParseTupleAndKeywords(args, kwargs, "ddddiidfOOOOO",
if (!PyArg_ParseTupleAndKeywords(args, kwargs, "ddddiidfOOOOOi",
keywords, &min_inclination, &max_inclination, &min_distance,
&max_distance, &prior_distance_power, &cosmology, &gmst, &sample_rate,
&epochs_obj, &snrs_obj, &responses_obj, &locations_obj, &horizons_obj))
&epochs_obj, &snrs_obj, &responses_obj, &locations_obj, &horizons_obj, &refine_level))
return NULL;
#pragma GCC diagnostic pop
......@@ -829,7 +830,7 @@ static PyObject *sky_map_toa_phoa_snr(
pixels = bayestar_sky_map_toa_phoa_snr(&len, &log_bci, &log_bsn,
min_inclination, max_inclination, min_distance, max_distance,
prior_distance_power, cosmology, gmst, nifos, nsamples, sample_rate,
epochs, snrs, responses, locations, horizons);
epochs, snrs, responses, locations, horizons, refine_level);
Py_END_ALLOW_THREADS
gsl_set_error_handler(old_handler);
......
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