diff --git a/.gitignore b/.gitignore index 0db678f8..009c30af 100644 --- a/.gitignore +++ b/.gitignore @@ -43,4 +43,5 @@ vlf_table.fits tests/.oda-token .ipynb_checkpoints -.venv \ No newline at end of file +cache +.venv diff --git a/oda_api/api.py b/oda_api/api.py index 8c1935d0..64ef5f35 100644 --- a/oda_api/api.py +++ b/oda_api/api.py @@ -2,6 +2,7 @@ from collections import OrderedDict from json.decoder import JSONDecodeError +from typing import Any, Callable, List, Tuple, Dict from astropy.table import Table from astropy.coordinates import Angle @@ -15,11 +16,11 @@ from .data_products import NumpyDataProduct, BinaryData, ApiCatalog, GWContoursDataProduct -from builtins import (bytes, str, open, super, range, - zip, round, input, int, pow, object, map, zip) __author__ = "Andrea Tramacere, Volodymyr Savchenko" +import hashlib +import gzip import warnings import requests import ast @@ -195,6 +196,9 @@ class DispatcherAPI: # allowing token discovery by default changes the user interface in some cases, # but in desirable way token_discovery_methods = None + use_local_cache = False + skip_parameter_check = False + raise_on_failure = True def __init__(self, instrument='mock', @@ -374,6 +378,12 @@ def selected_request_method(self): return self.preferred_request_method def request_to_json(self, verbose=False): + try: + return self.load_result() + except Exception as e: + logger.debug('unable to load result from %s: will need to compute', self.unique_response_json_fn) + + self.progress_logger.info( f'- waiting for remote response (since {time.strftime("%Y-%m-%d %H:%M:%S")}), please wait for {self.url}/{self.run_analysis_handle}') @@ -458,6 +468,9 @@ def request_to_json(self, verbose=False): self.returned_analysis_parameters = response_json['products'].get('analysis_parameters', None) + if response_json.get('query_status') in ['done', 'failed']: + self.save_result(response_json) + return response_json except json.decoder.JSONDecodeError as e: self.logger.error( @@ -662,7 +675,8 @@ def poll(self, verbose=None, silent=None): if self.is_complete: # TODO: something raising here does not help self.logger.debug("poll returing data: complete") - return DataCollection.from_response_json(self.response_json, self.instrument, self.product) + self.stored_result = DataCollection.from_response_json(self.response_json, self.instrument, self.product) + return self.stored_result def show_progress(self): full_report_dict_list = self.response_json['job_monitor'].get( @@ -741,10 +755,14 @@ def process_failure(self): self.response_json['exit_status']['message']) logger.error("have exception message: keys \"%s\"", exception_by_message.keys()) - raise exception_by_message.get(self.response_json['exit_status']['message'], RemoteException)( - message=self.response_json['exit_status']['message'], - debug_message=self.response_json['exit_status']['error_message'] - ) + + if self.raise_on_failure: + raise exception_by_message.get(self.response_json['exit_status']['message'], RemoteException)( + message=self.response_json['exit_status']['message'], + debug_message=self.response_json['exit_status']['error_message'] + ) + else: + self.exception_json = self.response_json['exit_status'] def failure_report(self, res_json): self.logger.error('query failed!') @@ -874,6 +892,42 @@ def report_last_request(self): self.logger.info( f"{C.GREY}last request completed in {self.last_request_t_complete - self.last_request_t0} seconds{C.NC}") + + def parameter_check(self, instrument, product, kwargs): + + res = requests.get("%s/api/par-names" % self.url, params=dict( + instrument=instrument, product_type=product), cookies=self.cookies) + + if res.status_code != 200: + warnings.warn( + 'parameter check not available on remote server, check carefully parameters name') + else: + _ignore_list = ['instrument', 'product_type', 'query_type', + 'off_line', 'query_status', 'verbose', 'session_id'] + validation_dict = copy.deepcopy(kwargs) + + for _i in _ignore_list: + del validation_dict[_i] + + valid_names = self._decode_res_json(res) + for n in validation_dict.keys(): + if n not in valid_names: + if self.strict_parameter_check: + raise UserError(f'the parameter: {n} is not among the valid ones: {valid_names}' + f'(you can set {self}.strict_parameter_check=False, but beware!') + else: + msg = '\n' + msg += '----------------------------------------------------------------------------\n' + msg += 'the parameter: %s ' % n + msg += ' is not among valid ones:' + msg += '\n' + msg += '%s' % valid_names + msg += '\n' + msg += 'this will throw an error in a future version \n' + msg += 'and might break the current request!\n ' + msg += '----------------------------------------------------------------------------\n' + warnings.warn(msg) + def get_list_terms_gallery(self, group: str = None, parent: str = None, @@ -1141,38 +1195,8 @@ def get_product(self, 'However the oda_api will perform a check of the list of valid parameters for your request.') del kwargs['dry_run'] - res = requests.get("%s/api/par-names" % self.url, params=dict( - instrument=instrument, product_type=product), cookies=self.cookies) - - if res.status_code != 200: - warnings.warn( - 'parameter check not available on remote server, check carefully parameters name') - else: - _ignore_list = ['instrument', 'product_type', 'query_type', - 'off_line', 'query_status', 'verbose', 'session_id'] - validation_dict = copy.deepcopy(kwargs) - - for _i in _ignore_list: - del validation_dict[_i] - - valid_names = self._decode_res_json(res) - for n in validation_dict.keys(): - if n not in valid_names: - if self.strict_parameter_check: - raise UserError(f'the parameter: {n} is not among the valid ones: {valid_names}' - f'(you can set {self}.strict_parameter_check=False, but beware!') - else: - msg = '\n' - msg += '----------------------------------------------------------------------------\n' - msg += 'the parameter: %s ' % n - msg += ' is not among valid ones:' - msg += '\n' - msg += '%s' % valid_names - msg += '\n' - msg += 'this will throw an error in a future version \n' - msg += 'and might break the current request!\n ' - msg += '----------------------------------------------------------------------------\n' - warnings.warn(msg) + if not self.skip_parameter_check: + self.parameter_check(instrument, product, kwargs) if kwargs.get('token', None) is None and self.token_discovery_methods is not None: discovered_token = oda_api.token.discover_token(self.token_discovery_methods) @@ -1202,12 +1226,10 @@ def get_product(self, d = DataCollection.from_response_json( res_json, instrument, product) - del (res) - return d @staticmethod - def set_api_code(query_dict, url="www.astro.unige.ch/mmoda/dispatch-data"): + def set_api_code(query_dict, url="www.astro.unige.ch/mmoda/dispatch-data") -> str: query_dict = OrderedDict(sorted(query_dict.items())) @@ -1247,14 +1269,45 @@ def __repr__(self): return f"[ {self.__class__.__name__}: {self.url} ]" -class DataCollection(object): + def save_result(self, response_json): + fn = self.unique_response_json_fn + + os.makedirs(os.path.dirname(fn), exist_ok=True) + + json.dump(response_json, gzip.open(fn, "wt")) + logger.info('saved result in %s', fn) + + + def load_result(self): + fn = self.unique_response_json_fn + logger.info('trying to load result from %s', fn) + + t0 = time.time() + r = json.load(gzip.open(fn, 'rb')) + + logger.info('\033[32mmanaged to load result\033[0m from %s in %.2f seconds', fn, time.time() - t0) + return r + + + @property + def unique_response_json_fn(self): + request_hash = hashlib.md5(self.set_api_code(self.parameters_dict).encode()).hexdigest()[:16] + + return f"cache/oda_api_data_collection_{request_hash}.json.gz" + + +class DataCollection(object): + + def __init__(self, data_list, add_meta_to_name=['src_name', 'product'], instrument=None, product=None): + self._p_list = [] + self._n_list = [] def __init__(self, data_list, add_meta_to_name=['src_name', 'product'], instrument=None, product=None, request_job_id=None): self._p_list = [] self._n_list = [] self.request_job_id = request_job_id + for ID, data in enumerate(data_list): - name = '' if hasattr(data, 'name'): name = data.name @@ -1296,10 +1349,42 @@ def as_list(self): meta_data = '' L.append({ - 'ID': ID, 'prod_name': prod_name, 'meta_data:': meta_data + 'ID': ID, + 'prod_name': prod_name, + 'metadata': meta_data, + 'meta_data:': meta_data # ??? }) - return L + return L + + + @property + def product_indexer(self): + return getattr(self, + '_product_indexer', + lambda p:(p['metadata']['src_name'], p['metadata']['product']) + ) + + + @product_indexer.setter + def product_indexer(self, v): + self._product_indexer = v + if len(self.keys()) != len(self.as_list()): + raise RuntimeError("duplicate index in metadata, this should be impossible, please check if product_indexer if good!") + + + def __getitem__(self, key): + return getattr(self, self.as_dict()[key]['prod_name']) + + + def as_dict(self): + return { self.product_indexer(product): product + for product in self.as_list() } + + + def keys(self): + return list(self.as_dict().keys()) + def _build_prod_name(self, prod, name, add_meta_to_name): @@ -1321,10 +1406,13 @@ def save_all_data(self, prenpend_name=None): file_name = file_name + '.fits' prod.write_fits_file(file_name) + def save(self, file_name): - pickle.dump(self, open(file_name, 'wb'), + pickle.dump(self, + gzip.open(file_name, 'wb'), protocol=pickle.HIGHEST_PROTOCOL) + def new_from_metadata(self, key, val): dc = None _l = [] @@ -1336,6 +1424,7 @@ def new_from_metadata(self, key, val): dc = DataCollection(_l) return dc + @classmethod def from_response_json(cls, res_json, instrument, product): @@ -1401,3 +1490,64 @@ def from_response_json(cls, res_json, instrument, product): p.meta_data = p.meta return d + + +class DispatcherAPICollection: + + def __init__(self, wait_between_poll_sequences_s=None, **kwargs) -> None: + self.wait_between_poll_sequences_s = wait_between_poll_sequences_s + self.constructor_kwargs = kwargs + + def get_product_list( + self, + parameter_dict_list: List[Dict[str, Any]], + **kwargs): + + self.client_list = [] + product_list = [] + + for parameter_dict in parameter_dict_list: + disp = DispatcherAPI(**self.constructor_kwargs) + + disp.use_local_cache = True + disp.wait = False + disp.skip_parameter_check = True + disp.raise_on_failure = False + + product_list.append(disp.get_product(**kwargs, **parameter_dict)) + + self.client_list.append(disp) + + logger.info('prepared %s clients, %s are done', + len(self.client_list), + len([c for c in self.client_list if c.is_complete]) + ) + + if self.wait_between_poll_sequences_s is None: + logger.info('not waiting for poll, please do not forget to come back for your results!') + else: + while True: + product_list = [] + + for client in self.client_list: + client.poll() + + logger.info('polled %s clients, %s are done', + len(self.client_list), + len([c for c in self.client_list if c.is_complete]), + ) + + product_list.append(client.poll()) + + if all([c.is_complete for c in self.client_list]): + logger.info('all done!') + break + else: + logger.info('will sleep %s s', + self.wait_between_poll_sequences_s + ) + time.sleep(self.wait_between_poll_sequences_s) + + return product_list + + diff --git a/oda_api/cli.py b/oda_api/cli.py index 48d8c5ea..169bc895 100644 --- a/oda_api/cli.py +++ b/oda_api/cli.py @@ -1,8 +1,6 @@ from datetime import datetime from email.policy import default import json -from attr import validate -from black import out import click import logging import time diff --git a/oda_api/data_products.py b/oda_api/data_products.py index 70682bca..03a71844 100644 --- a/oda_api/data_products.py +++ b/oda_api/data_products.py @@ -610,10 +610,52 @@ def decode(cls, encoded_obj: typing.Union[str, dict], from_json=False): -class ApiCatalog(object): +class ApiCatalog: - - def __init__(self,cat_dict,name='catalog'): + @classmethod + def from_list_of_dicts(cls, list_of_dicts): + { + "cat_column_descr": [ + ["meta_ID","f4"], + ["ra",">f4"], + ["dec",">f4"], + ["NEW_SOURCE",">i2"], + ["ISGRI_FLAG"," 1: + raise RuntimeError(f'object name {object_name} resolves to multiple different coordinates: {R}') + elif len(R) == 0: + raise RuntimeError(f'object name {object_name} can not be resolved with Simbad') + else: + source_coord = SkyCoord(R['RA'][0], R['DEC'][0], unit=(u.hourangle, u.deg)) # pylint: disable=no-member + return position_scw_list(source_coord, **kwargs) + + +def position_scw_list(source_coord, + radius_deg=5., + time_range="2001-01-01 .. 2029-01-01", + min_good_isgri=1000, + resultmax=int(1e6)): + + with astroquery.heasarc.Conf.server.set_temp('https://www.isdc.unige.ch/browse/w3query.pl'): + R = astroquery.heasarc.Heasarc().query_region( # pylint: disable=no-member + position=source_coord, + mission='integral_rev3_scw', resultmax=resultmax, radius=radius_deg*u.deg, cache=False, # pylint: disable=no-member + time=time_range, + fields='All', + good_isgri=f">{min_good_isgri}", + scw_type="POINTING" + + ) + + logger.debug('found %s SCWs', len(R)) + + return R \ No newline at end of file diff --git a/oda_api/tools/plot.py b/oda_api/tools/plot.py new file mode 100644 index 00000000..08fed50f --- /dev/null +++ b/oda_api/tools/plot.py @@ -0,0 +1,646 @@ +# mypy: ignore-errors +# pylint: skip-file +# pylint: disable-all + +from __future__ import absolute_import, division, print_function + +from builtins import (str, open, range, + zip, round, input, int, pow, object, zip) + + +__author__ = "Carlo Ferrigno" + +import json + +import numpy +from matplotlib import pylab as plt +from matplotlib.widgets import Slider, Button, RadioButtons +from matplotlib import cm + +import astropy.wcs as wcs +from astropy import table +from astropy import units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astroquery.simbad import Simbad +import copy + +import logging + +logger = logging.getLogger("oda_api.plot_tools") + +__all__ = ['OdaImage', 'OdaLightCurve'] + + +class OdaProduct(object): + + def __init__(self, data): + self.data = data + self.meta = None + self.logger = logger.getChild(self.__class__.__name__.lower()) + self.progress_logger = self.logger.getChild("progress") + + +class OdaImage(OdaProduct): + + def show(self, data=None, meta=None, header=None, sources=None, + levels=None, cmap=cm.gist_earth, + unit_ID=4, det_sigma=3): + + if levels is None: + levels = numpy.linspace(1, 10, 10) + + if data is None: + data = self.data.mosaic_image_0_mosaic.data_unit[unit_ID].data + + if meta is None: + self.meta = self.data.mosaic_image_0_mosaic.meta_data + + if header is None: + header = self.data.mosaic_image_0_mosaic.data_unit[unit_ID].header + + if sources is None: + sources = self.data.dispatcher_catalog_1.table + + fig = plt.figure(figsize=(8, 6)) + + j, i = plt.meshgrid(range(data.shape[0]), range(data.shape[1])) + w = wcs.WCS(header) + ra, dec = w.wcs_pix2world(numpy.column_stack([i.flatten(), j.flatten()]), 0).transpose() + ra = ra.reshape(i.shape) + dec = dec.reshape(j.shape) + + data = numpy.transpose(data) + data = numpy.ma.masked_equal(data, numpy.NaN) + + zero_crossing = False + + if numpy.abs(ra.max() - 360.0) < 0.1 and numpy.abs(ra.min()) < 0.1: + zero_crossing = True + ind_ra = ra > 180. + ra[ind_ra] -= 360. + ind_sort = numpy.argsort(ra, axis=-1) + ra = numpy.take_along_axis(ra, ind_sort, axis=-1) + data = numpy.take_along_axis(data, ind_sort, axis=-1) + + self.cs = plt.contourf(ra, dec, data, cmap=cmap, levels=levels, + extend="both", zorder=0) + self.cs.cmap.set_under('k') + self.cs.set_clim(numpy.min(levels), numpy.max(levels)) + + self.cb = plt.colorbar(self.cs) + + plt.xlim([ra.max(), ra.min()]) + plt.ylim([dec.min(), dec.max()]) + + if len(sources) > 0: + ras = numpy.array([x for x in sources['ra']]) + decs = numpy.array([x for x in sources['dec']]) + names = numpy.array([x for x in sources['src_names']]) + sigmas = numpy.array([x for x in sources['significance']]) + + # Defines relevant indexes for plotting regions + m_new = numpy.array(['NEW' in name for name in names]) + + # plot new sources as pink circles + try: + m = m_new & (sigmas > det_sigma) + ra_coord = ras[m] + dec_coord = decs[m] + new_names = names[m] + if zero_crossing: + ind_ra = ra_coord > 180. + try: + ra_coord[ind_ra] -= 360. + except: + pass + except: + ra_coord = [] + dec_coord = [] + new_names = [] + + plt.scatter(ra_coord, dec_coord, s=100, marker="o", facecolors='none', + edgecolors='pink', + lw=3, label="NEW any", zorder=5) + + for i in range(len(ra_coord)): + plt.text(ra_coord[i], + dec_coord[i] + 0.5, + new_names[i], color="pink", size=15) + + try: + m = ~m_new & (sigmas > det_sigma - 1) + ra_coord = ras[m] + dec_coord = decs[m] + cat_names = names[m] + if zero_crossing: + ind_ra = ra_coord > 180. + try: + ra_coord[ind_ra] -= 360. + except: + pass + except: + ra_coord = [] + dec_coord = [] + cat_names = [] + + plt.scatter(ra_coord, dec_coord, s=100, marker="o", facecolors='none', + edgecolors='magenta', lw=3, label="known", zorder=5) + + for i in range(len(ra_coord)): + plt.text(ra_coord[i], + dec_coord[i] + 0.5, + cat_names[i], color="magenta", size=15) + + plt.grid(color="grey", zorder=10) + + plt.xlabel("RA") + plt.ylabel("Dec") + + #Nice to have : slider + cmin = plt.axes([0.85, 0.05, 0.02, 0.4]) + cmax = plt.axes([0.85, 0.55, 0.02, 0.4]) + + data_min = data[numpy.isfinite(data)].min() + data_max = data[numpy.isfinite(data)].max() + + self.smin = Slider(cmin, 'Min', data_min, data_max, valinit=1., orientation='vertical') + self.smax = Slider(cmax, 'Max', data_min, data_max, valinit=10., orientation='vertical') + self.smin.on_changed(self.update) + self.smax.on_changed(self.update) + + plt.show() + + return fig + + def update(self, x): + if self.smin.val < self.smax.val: + self.cs.set_clim(self.smin.val, self.smax.val) + + + def write_fits(self, file_prefix=''): + self.data.mosaic_image_0_mosaic.write_fits_file(f'{file_prefix}mosaic.fits', overwrite=True) + + + def extract_catalog_from_image(self, include_new_sources=False, det_sigma=5, objects_of_interest=[], + flag=1, isgri_flag=2, update_catalog=False): + catalog_str = self.extract_catalog_string_from_image(include_new_sources, det_sigma, objects_of_interest, + flag, isgri_flag, update_catalog) + return json.loads(catalog_str) + + + def extract_catalog_string_from_image(self, include_new_sources=False, det_sigma=5, + objects_of_interest=None, + flag=1, isgri_flag=2, update_catalog=True) -> str: + """ + Example: objects_of_interest=['Her X-1'] + objects_of_interest=[('Her X-1', Simbad.query )] + objects_of_interest=[('Her X-1', Skycoord )] + objects_of_interest=[ Skycoord(....) ] + """ + + if objects_of_interest is None: + objects_of_interest = [] + + image = self.data + + if image.dispatcher_catalog_1.table is None: + self.logger.warning("No sources in the catalog") + if objects_of_interest != []: + return OdaImage.add_objects_of_interest(None, objects_of_interest, + flag, isgri_flag) + else: + return 'none' + + sources = image.dispatcher_catalog_1.table[image.dispatcher_catalog_1.table['significance'] >= det_sigma] + + if len(sources) == 0: + self.logger.warning('No sources in the catalog with det_sigma > %.1f' % det_sigma) + if objects_of_interest != []: + return self.add_objects_of_interest(None, objects_of_interest, + flag, isgri_flag) + else: + return 'none' + + if not include_new_sources: + ind = [not 'NEW' in ss for ss in sources['src_names']] + clean_sources = sources[ind] + self.logger.debug(ind) + self.logger.debug(sources) + self.logger.debug(clean_sources) + else: + clean_sources = sources + + unique_sources = self.add_objects_of_interest(clean_sources, objects_of_interest, + flag, isgri_flag) + + copied_image = copy.deepcopy(image) + copied_image.dispatcher_catalog_1.table = unique_sources + + if update_catalog: + image.dispatcher_catalog_1.table = unique_sources + + return copied_image.dispatcher_catalog_1.get_api_dictionary() + + @staticmethod + def make_one_source_catalog_string(name, ra, dec, isgri_flag, flag): + out_str_templ ='{"cat_frame": "fk5", "cat_coord_units": "deg", "cat_column_list": [[1], ["%s"], [0.0], [%f], [%f], [-32768], [%d], [%d], [0.001]], "cat_column_names": ["meta_ID", "src_names", "significance", "ra", "dec", "NEW_SOURCE", "ISGRI_FLAG", "FLAG", "ERR_RAD"], "cat_column_descr": [["meta_ID", " 0: + self.logger.info('Found ' + ooi + ' in catalog') + clean_sources['FLAG'][ind] = flag + if 'ISGRI_FLAG' in clean_sources.keys(): + clean_sources['ISGRI_FLAG'][ind] = isgri_flag + if 'JEMX_FLAG' in clean_sources.keys(): + clean_sources['JEMX_FLAG'][ind] = isgri_flag + else: + self.logger.info('Adding ' + ooi + ' to catalog') + try: + self.logger.debug('Flux is present') + clean_sources.add_row((0, ooi, 0, ra, dec, 0, isgri_flag, flag, 1e-3, 0, 0)) + except: + self.logger.debug('Flux is NOT present') + clean_sources.add_row((0, ooi, 0, ra, dec, 0, isgri_flag, flag, 1e-3)) + + unique_sources = table.unique(clean_sources, keys=['src_names']) + + return unique_sources + else: + return self.make_one_source_catalog_string(ooi, ra, dec, isgri_flag, flag) + + +class OdaLightCurve(OdaProduct): + + def get_lc(self, source_name, systematic_fraction=0): + + combined_lc = self.data + # In LC name has no "-" nor "+" ?????? + patched_source_name = source_name.replace('-', ' ').replace('+', ' ') + + hdu = None + for j, dd in enumerate(combined_lc._p_list): + self.logger.debug(dd.meta_data['src_name']) + if dd.meta_data['src_name'] == source_name or dd.meta_data['src_name'] == patched_source_name: + for ii, du in enumerate(dd.data_unit): + if 'LC' in du.name: + hdu = du.to_fits_hdu() + + if hdu is None: + self.logger.info('Source ' + source_name + ' not found in the light curves') + return None, None, None, None, None, None + + x = hdu.data['TIME'] + y = hdu.data['RATE'] + dy = hdu.data['ERROR'] + self.logger.debug("Original length of light curve %d" % len(x)) + ind = numpy.argsort(x) + x = x[ind] + y = y[ind] + dy = dy[ind] + dy = numpy.sqrt(dy ** 2 + (y * systematic_fraction) ** 2) + ind = numpy.logical_and(numpy.isfinite(y), numpy.isfinite(dy)) + ind = numpy.logical_and(ind, dy > 0) + self.logger.debug("Final length of light curve %d " % numpy.sum(ind)) + + try: + e_min = hdu.header['E_MIN'] + except: + e_min = 0 + + try: + e_max = hdu.header['E_MAX'] + except: + e_max = 0 + + #This could only be valid for ISGRI + try: + dt_lc = hdu.data['XAX_E'] + self.logger.debug('Get time bin directly from light curve') + except: + timedel = hdu.header['TIMEDEL'] + timepix = hdu.header['TIMEPIXR'] + t_lc = hdu.data['TIME'] + (0.5 - timepix) * timedel + dt_lc = t_lc.copy() * 0.0 + timedel / 2 + for i in range(len(t_lc) - 1): + dt_lc[i + 1] = min(timedel / 2, t_lc[i + 1] - t_lc[i] - dt_lc[i]) + self.logger.debug('Computed time bin from TIMEDEL') + + return x[ind], dt_lc[ind], y[ind], dy[ind], e_min, e_max + + def show(self, in_source_name='', systematic_fraction=0, ng_sig_limit=0, find_excesses=False): + #if ng_sig_limit <1 does not plot range + combined_lc = self.data + from scipy import stats + + if in_source_name == '': + source_names = [dd.meta_data['src_name'] for dd in combined_lc._p_list] + else: + source_names = [in_source_name] + + for source_name in source_names: + x, dx, y, dy, e_min, e_max = self.get_lc(source_name, systematic_fraction) + if x is None: + return + + meany = numpy.sum(y / dy ** 2) / numpy.sum(1. / dy ** 2) + err_mean = numpy.sum(1 / dy ** 2) + + std_dev = numpy.std(y) + + fig = plt.figure() + _ = plt.errorbar(x, y, xerr=dx, yerr=dy, marker='o', capsize=0, linestyle='', label='Lightcurve') + _ = plt.axhline(meany, color='green', linewidth=3) + _ = plt.xlabel('Time [IJD]') + if e_min == 0 or e_max ==0: + _ = plt.ylabel('Rate') + else: + _ = plt.ylabel('Rate %.1f-%.1f keV' % (e_min, e_max)) + + if ng_sig_limit >= 1: + ndof = len(y) - 1 + prob_limit = stats.norm().sf(ng_sig_limit) + chi2_limit = stats.chi2(ndof).isf(prob_limit) + band_width = numpy.sqrt(chi2_limit / err_mean) + + _ = plt.axhspan(meany - band_width, meany + band_width, color='green', alpha=0.3, + label=f'{ng_sig_limit} $\sigma_m$, {100 * systematic_fraction}% syst') + + _ = plt.axhspan(meany - std_dev*ng_sig_limit, meany + std_dev*ng_sig_limit, + color='cyan', alpha=0.3, + label=f'{ng_sig_limit} $\sigma_d$, {100 * systematic_fraction}% syst') + + _ = plt.legend() + + plot_title = source_name + _ = plt.title(plot_title) + if find_excesses: + ind = (y - band_width)/dy > ng_sig_limit + if numpy.sum(ind) > 0: + _ = plt.plot(x[ind], y[ind], marker='x', color='red', linestyle='', markersize=10) + self.logger.info('We found positive excesses on the lightcurve at times') + good_ind = numpy.where(ind) + #print(good_ind[0][0:-1], good_ind[0][1:]) + old_time = -1 + if len(good_ind[0]) == 1: + self.logger.info('%f' % (x[good_ind[0][0]])) + else: + for i,j in zip(good_ind[0][0:-1], good_ind[0][1:]): + #print(i,j) + if j-i > 2: + if x[i] != old_time : + self.logger.info('%f' % x[i]) + _ = OdaLightCurve.plot_zoom(x,y,dy,i) + self.logger.info('%f' % (x[j])) + _ = OdaLightCurve.plot_zoom(x, y, dy, j) + # else: + # self.logger.debug('%f' % ((x[i]+x[j])/2)) + + old_time = x[j] + + return fig + + + @staticmethod + def plot_zoom(x, y, dy, i, n_before=5, n_after=15, save_plot=True, name_base='burst_at_'): + fig = plt.figure() + _ = plt.errorbar(x[i-n_before:i+n_after], y[i-n_before:i+n_after], yerr=dy[i-n_before:i+n_after], + marker='o', capsize=0, linestyle='', label='Lightcurve') + _ = plt.xlabel('Time [IJD]') + _ = plt.ylabel('Rate') + if save_plot: + _ = plt.savefig(name_base+'%d.png' % i) + return fig + + + def write_fits(self, source_name, file_suffix='', output_dir='.'): + # In LC name has no "-" nor "+" ?????? + lc = self.data + patched_source_name = source_name.replace('-', ' ').replace('+', ' ') + lcprod = [l for l in lc._p_list if l.meta_data['src_name'] == source_name or \ + l.meta_data['src_name'] == patched_source_name] + + if (len(lcprod) < 1): + self.logger.warning("source %s not found in light curve products" % source_name) + return "none", 0, 0, 0 + + if (len(lcprod) > 1): + self.logger.warning( + "source %s is found more than once light curve products, writing only the first one" % source_name) + + instrument = lcprod[0].data_unit[1].header['INSTRUME'] + if instrument == 'IBIS': + ind_extension = 1 + else: + ind_extension = 2 + + lc_fn = output_dir + "/%s_lc_%s%s.fits" % (instrument, source_name.replace(' ', '_'), file_suffix) + hdu = lcprod[0].data_unit[ind_extension].to_fits_hdu() + timedel = hdu.header['TIMEDEL'] + timepixr = hdu.header['TIMEPIXR'] + + dt = timedel * timepixr + + hdu.header['TSTART'] = hdu.data['TIME'][0] - dt + hdu.header['TSTOP'] = hdu.data['TIME'][-1] + dt + hdu.header['TFIRST'] = hdu.data['TIME'][0] - dt + hdu.header['TLAST'] = hdu.data['TIME'][-1] + dt + hdu.header['TELAPSE'] = hdu.header['TLAST'] - hdu.header['TFIRST'] + + ontime=0 + for x in hdu.data['FRACEXP']: + ontime += x * timedel + + hdu.header['ONTIME'] = ontime + + fits.writeto(lc_fn, hdu.data, header=hdu.header, overwrite=True) + + mjdref = float(hdu.header['MJDREF']) + tstart = float(hdu.header['TSTART']) + mjdref + tstop = float(hdu.header['TSTOP']) + mjdref + try: + exposure = float(hdu.header['EXPOSURE']) + except: + exposure = -1 + + return lc_fn, tstart, tstop, exposure + + +class OdaSpectrum(OdaProduct): + + def show_spectral_products(self): + + summed_data = self.data + + for dd, nn in zip(summed_data._p_list, summed_data._n_list): + self.logger.debug(nn) + dd.show_meta() + # for kk in dd.meta_data.items(): + if 'spectrum' in dd.meta_data['product']: + self.logger.debug(dd.data_unit[1].header['EXPOSURE']) + dd.show() + + def get_spectrum_products(self, in_source_name='none'): + if in_source_name == 'none': + return None + + specprod = [l for l in self.data._p_list if l.meta_data['src_name'] == in_source_name] + + if (len(specprod) < 1): + self.logger.warning("source %s not found in spectral products" % in_source_name) + return None + + return specprod + + def show(self, in_source_name='', systematic_fraction=0, xlim=[]): + + if in_source_name == '': + self.show_spectral_products() + return + + specprod = self.get_spectrum_products(in_source_name) + if specprod is None: + return + + spec = specprod[0].data_unit[1].to_fits_hdu() + for hh in specprod[2].data_unit: + if hh.to_fits_hdu().header['EXTNAME'] == 'EBOUNDS': + ebounds = hh.to_fits_hdu() + + x = (ebounds.data['E_MAX'] + ebounds.data['E_MIN'])/2. + dx = (ebounds.data['E_MAX'] - ebounds.data['E_MIN']) / 2. + y = spec.data['RATE'] + dy = numpy.sqrt(spec.data['STAT_ERR']**2 + spec.data['SYS_ERR']**2 + (y*systematic_fraction)**2) + + fig = plt.figure() + _ = plt.errorbar(x, y, xerr=dx, yerr=dy, marker='o', capsize=0, linestyle='', label='spectrum') + + _ = plt.xlabel('Energy [keV]') + _ = plt.xscale('log') + _ = plt.yscale('log') + _ = plt.ylabel('Rate') + _ = plt.title(in_source_name) + if len(xlim) == 2: + _ = plt.xlim(xlim) + + return fig + + def write_fits(self, source_name='', file_suffix='', grouping=[0, 0, 0], systematic_fraction=0, + output_dir='.'): + """ + Grouping argument is [minimum_energy, maximum_energy, number_of_bins] + number of bins > 0, linear grouping + number_of_bins < 0, logarithmic binning + """ + + if source_name == '': + self.show_spectral_products() + self.logger.warning('PLease specify a source to save the spectral products') + return "none", 0, 0, 0 + + specprod = self.get_spectrum_products(source_name) + if specprod is None: + return "none", 0, 0, 0 + + instrument = specprod[0].data_unit[1].header['INSTRUME'] + + out_name = source_name.replace(' ', '_').replace('+', 'p') + spec_fn = output_dir + "/%s_spectrum_%s%s.fits" % (instrument, out_name, file_suffix) + arf_fn = output_dir + "/%s_arf_%s%s.fits" % (instrument, out_name, file_suffix) + rmf_fn = output_dir + "/%s_rmf_%s%s.fits" % (instrument, out_name, file_suffix) + + self.logger.info("Saving spectrum %s with rmf %s and arf %s" % (spec_fn, rmf_fn, arf_fn)) + + specprod[0].write_fits_file(spec_fn) + specprod[1].write_fits_file(arf_fn) + specprod[2].write_fits_file(rmf_fn) + + ff = fits.open(spec_fn, mode='update') + + ff[1].header['RESPFILE'] = rmf_fn + ff[1].header['ANCRFILE'] = arf_fn + mjdref = ff[1].header['MJDREF'] + tstart = float(ff[1].header['TSTART']) + mjdref + tstop = float(ff[1].header['TSTOP']) + mjdref + exposure = ff[1].header['EXPOSURE'] + ff[1].data['SYS_ERR'] = numpy.zeros(len(ff[1].data['SYS_ERR'])) + systematic_fraction + ind = numpy.isfinite(ff[1].data['RATE']) + ff[1].data['QUALITY'][ind] = 0 + + if numpy.sum(grouping) != 0: + + if grouping[1] <= grouping[0] or grouping[2] == 0: + raise RuntimeError('Wrong grouping arguments') + + ff_rmf = fits.open(rmf_fn) + + e_min = ff_rmf['EBOUNDS'].data['E_MIN'] + e_max = ff_rmf['EBOUNDS'].data['E_MAX'] + + ff_rmf.close() + + ind1 = numpy.argmin(numpy.abs(e_min - grouping[0])) + ind2 = numpy.argmin(numpy.abs(e_max - grouping[1])) + + n_bins = numpy.abs(grouping[2]) + + ff[1].data['GROUPING'][0:ind1] = 0 + ff[1].data['GROUPING'][ind2:] = 0 + + ff[1].data['QUALITY'][0:ind1] = 1 + ff[1].data['QUALITY'][ind2:] = 1 + + if grouping[2] > 0: + step = int((ind2 - ind1 + 1) / n_bins) + self.logger.info('Linear grouping with step %d' % step) + for i in range(1, step): + j = range(ind1 + i, ind2, step) + ff[1].data['GROUPING'][j] = -1 + else: + ff[1].data['GROUPING'][ind1:ind2] = -1 + e_step = (e_max[ind2] / e_min[ind1]) ** (1.0 / n_bins) + self.logger.info('Geometric grouping with step %.3f' % e_step) + loc_e = e_min[ind1] + while (loc_e < e_max[ind2]): + ind_loc_e = numpy.argmin(numpy.abs(e_min - loc_e)) + ff[1].data['GROUPING'][ind_loc_e] = 1 + loc_e *= e_step + + ff.flush() + ff.close() + + return spec_fn, tstart, tstop, exposure + diff --git a/tests/test_basic.py b/tests/test_basic.py index dacd5aac..36550af2 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -376,4 +376,107 @@ def get_exception(*args, **kwargs): else: raise RuntimeError() - requests.get = requests._get \ No newline at end of file + requests.get = requests._get + + +def test_storing(dispatcher_api): + disp = dispatcher_api + + token_payload = { + **default_token_payload, + "roles": ["general", "integral-private-qla"], + } + encoded_token = jwt.encode(token_payload, secret_key, algorithm='HS256') + + dc = disp.get_product( + product_type="Dummy", + instrument="empty", + product="numerical", + token=encoded_token + ) + + disp.save_result() + + assert disp.load_result().as_list() == dc.as_list() + + +def test_multiple_requests(dispatcher_live_fixture): + import oda_api.api + + dac = oda_api.api.DispatcherAPICollection(url=dispatcher_live_fixture, + wait_between_poll_sequences_s=None) + + dac.get_product_list([ + dict(product="dummy", + instrument="empty", + RA=i, + p_list=['111100110010.001']) + for i in [1, 2] + ], + ) + + +@pytest.mark.isgri() +def test_indexing(dispatcher_api): + disp = dispatcher_api + + token_payload = { + **default_token_payload, + "roles": ["general", "integral-private-qla"], + } + encoded_token = jwt.encode(token_payload, secret_key, algorithm='HS256') + + dc = disp.get_product( + product_type="Dummy", + instrument="empty", + product="numerical", + token=encoded_token + ) + + assert isinstance(dc, oda_api.api.DataCollection) + + assert dc.as_list() == [{'ID': 0, 'prod_name': 'numerical_0', 'metadata': {}, 'meta_data:': {}}] + + dc.product_indexer = lambda p: p['prod_name'] + + assert dc.as_dict() == {'numerical_0': {'ID': 0, 'prod_name': 'numerical_0', 'metadata': {}, 'meta_data:': {}}} + + assert dc.keys() == ['numerical_0'] + + assert dc['numerical_0'] + + +def test_local_cache_unit(): + from oda_api.localcache import cached, call_to_fn + from collections import OrderedDict + + fn = 'cache-test.txt' + with open(fn, 'wt') as f: + f.write('test1') + + @cached + def testme(x, a): + return open(fn).read() + str(x) + str(a) + + cache_fn = call_to_fn(testme, 1, a='b') + + if os.path.exists(cache_fn): + os.remove(cache_fn) + + assert testme(1, a='b') == "test11b" + + with open(fn, 'wt') as f: + f.write('test2') + + # did not change since it's cached + assert testme(1, a='b') == "test11b" + + os.remove(call_to_fn(testme, 1, a='b')) + + assert testme(1, a='b') == "test21b" + assert testme(2, a='b') == "test22b" + + + +def test_local_cache(dispatcher_api): + pass \ No newline at end of file diff --git a/tests/test_tools_integral.py b/tests/test_tools_integral.py new file mode 100644 index 00000000..8dc99fed --- /dev/null +++ b/tests/test_tools_integral.py @@ -0,0 +1,11 @@ + +def test_scwlist(): + from astropy.coordinates import SkyCoord + from oda_api.tools import integral + R = integral.position_scw_list(SkyCoord(83, 22, unit='deg'), resultmax=10) + assert len(R) == 10 + +def test_object_scwlist(): + from oda_api.tools import integral + R = integral.object_scw_list('Crab', resultmax=10) + assert len(R) == 10 \ No newline at end of file