From 26f848d274f6eced3fd06d1a73d144607dec30d4 Mon Sep 17 00:00:00 2001 From: tyrin Date: Mon, 19 May 2025 15:09:04 +0300 Subject: [PATCH] code formatting --- management/commands/load_survey.py | 83 ++--- management/commands/set_contaminated.py | 216 +++++++------ models.py | 66 ++-- serializers.py | 5 - urls.py | 8 +- views.py | 412 +++++++++++------------- 6 files changed, 373 insertions(+), 417 deletions(-) diff --git a/management/commands/load_survey.py b/management/commands/load_survey.py index 553160b..f792ae9 100644 --- a/management/commands/load_survey.py +++ b/management/commands/load_survey.py @@ -1,21 +1,21 @@ # uplim/management/commands/load_survey.py import numpy as np + from astropy.io import fits +from itertools import islice +from datetime import datetime from django.core.management.base import BaseCommand from django.db import transaction from uplim.models import Pixel from django.db.models import Max -from itertools import islice - -from datetime import datetime - # DEFINE BATCH SIZE AND BATCH # ************************************************************** -#BATCH_SIZE = 1000000 +# BATCH_SIZE = 1000000 + def batch(iterable, size): """ @@ -29,8 +29,6 @@ def batch(iterable, size): yield chunk - - class Command(BaseCommand): help = "Process FITS files and store the data in the database" @@ -39,40 +37,33 @@ class Command(BaseCommand): def add_arguments(self, parser): parser.add_argument( - '--counts', - type=str, - required=True, - help='Path of the counts file' + "--counts", type=str, required=True, help="Path of the counts file" ) parser.add_argument( - '--exposure', - type=str, - required=True, - help='Path of the exposure file' + "--exposure", type=str, required=True, help="Path of the exposure file" ) parser.add_argument( - '--survey_number', - type=int, + "--survey_number", + type=int, required=True, - help='Integer ID of the survey being read' + help="Integer ID of the survey being read", ) parser.add_argument( - '--batch_size', - type=int, + "--batch_size", + type=int, default=1000, - help='Integer number of pixels to be inserted into the database at once' + help="Integer number of pixels to be inserted into the database at once", ) - def handle(self, *args, **options): - + # GET FILENAMES FROM ARGUMENTS # ************************************************************** - - counts_file = options['counts'] - exposure_file = options['exposure'] - survey_number = options['survey_number'] - BATCH_SIZE = options['batch_size'] + + counts_file = options["counts"] + exposure_file = options["exposure"] + survey_number = options["survey_number"] + BATCH_SIZE = options["batch_size"] self.stdout.write(f"\nCounts file:\t{counts_file}") self.stdout.write(f"Exposure file:\t{exposure_file}") @@ -87,7 +78,6 @@ class Command(BaseCommand): counts_data = counts_map.ravel() - with fits.open(exposure_file) as hdul: column_name = "T" @@ -100,49 +90,48 @@ class Command(BaseCommand): self.stdout.write(f"\nCounts Data Shape:\t{counts_data.shape}") self.stdout.write(f"Exposure Data Shape:\t{exposure_data.shape}") - + total_pixels = counts_data.shape[0] self.stdout.write(f"\nTotal pixels to insert:\t{total_pixels}") - assert counts_data.shape == exposure_data.shape, "Counts and exposure maps must have the same shape" + assert ( + counts_data.shape == exposure_data.shape + ), "Counts and exposure maps must have the same shape" # CREATE THE SURVEY IF IT DOES NOT EXIST # ************************************************************** - + # with transaction.atomic(): - + # survey,created = Survey.objects.get_or_create(number=survey_number) - + # if created: # self.stdout.write(f"Created a new survey instance with number: {survey.number}") # else: # self.stdout.write(f"Using existing survey instance with the number: {survey.number}") - + # FETCH THE LAST PROCESSED HPID AND CONTINUE FROM IT # ************************************************************** - + last_hpid = ( - Pixel.objects - .filter(survey=survey_number) - .aggregate(max_hpid=Max('hpid'))['max_hpid'] + Pixel.objects.filter(survey=survey_number).aggregate(max_hpid=Max("hpid"))[ + "max_hpid" + ] or -1 ) start_index = last_hpid + 1 - - - + pixel_generator = ( Pixel( hpid=i, counts=int(count), exposure=float(exposure), - survey=survey_number + survey=survey_number, ) for i, (count, exposure) in enumerate(zip(counts_data, exposure_data)) - if i >= start_index + if i >= start_index ) - total_inserted = start_index # Process in batches for pixel_batch in batch(pixel_generator, BATCH_SIZE): @@ -151,8 +140,6 @@ class Command(BaseCommand): total_inserted += len(pixel_batch) percentage = total_inserted / total_pixels * 100 timestamp = datetime.now().strftime("%H:%M:%S") - self.stdout.write( - f"[{timestamp}] {percentage:.2f}% inserted" - ) + self.stdout.write(f"[{timestamp}] {percentage:.2f}% inserted") self.stdout.write(f"Inserted a total of {total_inserted} pixels.") diff --git a/management/commands/set_contaminated.py b/management/commands/set_contaminated.py index 8f32a94..539b12f 100644 --- a/management/commands/set_contaminated.py +++ b/management/commands/set_contaminated.py @@ -7,19 +7,19 @@ from django.core.management.base import BaseCommand from django.db import transaction +from uplim.models import Pixel, CatalogSource import pandas as pd import healpy as hp import numpy as np from astropy.coordinates import SkyCoord -from uplim.models import Pixel, CatalogSource - from itertools import islice - from datetime import datetime -BATCH_SIZE=900 + +BATCH_SIZE = 900 + def batch(iterable, size): iterable = iter(iterable) @@ -28,157 +28,165 @@ def batch(iterable, size): if not chunk: break yield chunk - + class Command(BaseCommand): help = "Set the 'contaminated' flag for all pixels based on the fluxes in the provided catalog." - + # COMMAND LINE ARGUMENTS # ********************** - + def add_arguments(self, parser): - + parser.add_argument( - '--catalog', - type=str, - required=False, - help='Path to the catalog.dat file' + "--catalog", type=str, required=False, help="Path to the catalog.dat file" ) - + # parser.add_argument( # '--survey', # type=int, # required=False, # help='integer number of the survey to set the flag for' # ) - + parser.add_argument( - '--reset', - action='store_true', + "--reset", + action="store_true", default=False, - help='Reset the contamination flag across all pixels back to False.' + help="Reset the contamination flag across all pixels back to False.", ) - + def handle(self, *args, **options): - + # RESET BEHAVIOR: SET CONTAMINATION FLAG TO FALSE FOR ALL PIXELS # ************************************************************** - - if options['reset']: - + + if options["reset"]: + self.stdout.write("Resetting the contamination flag...") - - Pixel.objects.update(contaminated = False) - + + Pixel.objects.update(contaminated=False) + self.stdout.write("Done") return - - if not options['catalog']: + + if not options["catalog"]: self.stdout.write("No catalog file provided, exiting") return - - catalog_file = options['catalog'] - + + catalog_file = options["catalog"] + self.stdout.write(f"Catalog file:\t{catalog_file}") - + # READ THE CATALOG FILE USING PANDAS READ_FWF # ******************************************* - - # Define column positions based on the byte ranges + + # Define column positions based on the byte ranges colspecs = [ - (0, 4), # SrcID (1-4) - (5, 26), # Name (6-26) + (0, 4), # SrcID (1-4) + (5, 26), # Name (6-26) (27, 37), # RAdeg (28-37) (38, 48), # DEdeg (39-48) (49, 55), # ePos (50-55) (56, 63), # Signi (57-63) (64, 76), # Flux (65-76) (77, 89), # e_Flux (78-89) - (90, 118), # CName (91-118) - (119, 120),# NewXray (120) - (121, 134) # Type (122-134) + (90, 118), # CName (91-118) + (119, 120), # NewXray (120) + (121, 134), # Type (122-134) ] # Define column names colnames = [ - "SrcID", "Name", "RAdeg", "DEdeg", "ePos", "Signi", "Flux", - "e_Flux", "CName", "NewXray", "Type" + "SrcID", + "Name", + "RAdeg", + "DEdeg", + "ePos", + "Signi", + "Flux", + "e_Flux", + "CName", + "NewXray", + "Type", ] # Read the file using the fixed-width format catalog = pd.read_fwf(catalog_file, colspecs=colspecs, names=colnames) - - for col in ['Name', 'CName', 'Type']: - catalog[col] = catalog[col].fillna('') - + + for col in ["Name", "CName", "Type"]: + catalog[col] = catalog[col].fillna("") + self.stdout.write(str(catalog.head())) - + # LOAD THE CATALOG INTO THE DATABASE # ********************************** - - existing_srcids = set( - CatalogSource.objects.values_list('srcid', flat=True) - ) - + + existing_srcids = set(CatalogSource.objects.values_list("srcid", flat=True)) + to_create = [] - + for _, row in catalog.iterrows(): - - srcid = int(row['SrcID']) + + srcid = int(row["SrcID"]) if srcid in existing_srcids: continue to_create.append( CatalogSource( - srcid = srcid, - name = row['Name'].strip(), - ra_deg = float(row['RAdeg']), - dec_deg = float(row['DEdeg']), - pos_error = float(row['ePos']), - significance = float(row['Signi']), - flux = float(row['Flux']), - flux_error = float(row['e_Flux']), - catalog_name = row['CName'].strip(), - new_xray = bool(int(row['NewXray'])), - source_type = row['Type'].strip() + srcid=srcid, + name=row["Name"].strip(), + ra_deg=float(row["RAdeg"]), + dec_deg=float(row["DEdeg"]), + pos_error=float(row["ePos"]), + significance=float(row["Signi"]), + flux=float(row["Flux"]), + flux_error=float(row["e_Flux"]), + catalog_name=row["CName"].strip(), + new_xray=bool(int(row["NewXray"])), + source_type=row["Type"].strip(), ) ) - + if to_create: - self.stdout.write(f'Inserting {len(to_create)} new catalog rows.') + self.stdout.write(f"Inserting {len(to_create)} new catalog rows.") for chunk in batch(to_create, BATCH_SIZE): CatalogSource.objects.bulk_create(chunk, ignore_conflicts=True) - self.stdout.write('Catalog update complete.') + self.stdout.write("Catalog update complete.") else: - self.stdout.write('All catalog rows already exist in the database.') - + self.stdout.write("All catalog rows already exist in the database.") + # hard coded nside and flux-radius mapping # maybe change that - + nside = 4096 npix = hp.nside2npix(nside) - - flux_bins = [0, 125, 250, 2000, 20000, np.inf] # define bin edges - mask_radii_deg = [ 0.06, 0.15, 0.5, 0.9, 2.5 ] # corresponding mask radii in degrees - + + flux_bins = [0, 125, 250, 2000, 20000, np.inf] # define bin edges + mask_radii_deg = [ + 0.06, + 0.15, + 0.5, + 0.9, + 2.5, + ] # corresponding mask radii in degrees + # Convert mask radii from degrees to radians (required by query_disc) mask_radii = [np.radians(r) for r in mask_radii_deg] # Use pandas.cut to assign each source a bin index (0, 1, or 2) - catalog['flux_bin'] = pd.cut(catalog['Flux'], bins=flux_bins, labels=False) + catalog["flux_bin"] = pd.cut(catalog["Flux"], bins=flux_bins, labels=False) # manually add and change some sources manual_additions = pd.DataFrame( [ - {'RAdeg' : 279.9804336, 'DEdeg' : 5.0669542, 'flux_bin' : 3}, - {'RAdeg' : 266.5173685, 'DEdeg' : -29.1252321, 'flux_bin' : 3}, + {"RAdeg": 279.9804336, "DEdeg": 5.0669542, "flux_bin": 3}, + {"RAdeg": 266.5173685, "DEdeg": -29.1252321, "flux_bin": 3}, ] ) catalog = pd.concat([catalog, manual_additions], ignore_index=True) - catalog.loc[catalog['SrcID'] == 1101, 'flux_bin'] = 2 - - + catalog.loc[catalog["SrcID"] == 1101, "flux_bin"] = 2 mask_array = np.ones(npix, dtype=bool) @@ -188,62 +196,58 @@ class Command(BaseCommand): # process each source in the catalog for _, row in catalog.iterrows(): - - ra = row['RAdeg'] - dec = row['DEdeg'] - - src_coord = SkyCoord( - ra, dec, unit = 'deg', frame = 'icrs' - ) + + ra = row["RAdeg"] + dec = row["DEdeg"] + + src_coord = SkyCoord(ra, dec, unit="deg", frame="icrs") gal = src_coord.galactic - + ra, dec = gal.l.deg, gal.b.deg - - flux_bin = row['flux_bin'] # 0, 1, or 2 + + flux_bin = row["flux_bin"] # 0, 1, or 2 # Get the corresponding mask radius (in radians) for this flux bin radius = mask_radii[flux_bin] - + # Convert (ra, dec) to HEALPix spherical coordinates theta = np.radians(90.0 - dec) phi = np.radians(ra) vec = hp.ang2vec(theta, phi) - + # Query all pixels within the given radius # 'inclusive=True' makes sure pixels on the edge are included pix_indices = hp.query_disc(nside, vec, radius, inclusive=True) - + # Mark these pixels as bad (False) in our mask mask_array[pix_indices] = False # Add the pixel indices to our set of masked pixels masked_pixels_set.update(pix_indices) - - # Convert the set of masked pixels to a sorted list. masked_pixels_list = sorted(list(masked_pixels_set)) # print("Number of masked pixels:", len(masked_pixels_list)) self.stdout.write("\nList ready, updating the database...") - - - + if not masked_pixels_list: self.stdout.write("No pixels marked as contaminated, exiting.") return - + total = len(masked_pixels_list) updated = 0 - self.stdout.write(f'\nUpdating contaminated flag in batches of {BATCH_SIZE}') - + self.stdout.write(f"\nUpdating contaminated flag in batches of {BATCH_SIZE}") + for chunk in batch(masked_pixels_list, BATCH_SIZE): with transaction.atomic(): Pixel.objects.filter(hpid__in=chunk).update(contaminated=True) - + updated += len(chunk) percentage = updated / total * 100 - + timestamp = datetime.now().strftime("%H:%M:%S") - self.stdout.write(f'[{timestamp}] {updated}/{total} ({percentage:.1f}%) updated') - - self.stdout.write(f'\n Marked {updated} pixels as contaminated.') \ No newline at end of file + self.stdout.write( + f"[{timestamp}] {updated}/{total} ({percentage:.1f}%) updated" + ) + + self.stdout.write(f"\n Marked {updated} pixels as contaminated.") diff --git a/models.py b/models.py index 2c19b02..6b193c7 100644 --- a/models.py +++ b/models.py @@ -4,51 +4,47 @@ from django.db import models from django.db.models import UniqueConstraint - - class Pixel(models.Model): - - #id = models.AutoField(primary_key=True) # ~200 million pixels for a 4096 survey + + # id = models.AutoField(primary_key=True) # ~200 million pixels for a 4096 survey # no need to set explicitly # WILL ONLY HOLD 10 SURVEYS AS AN AUTOFIELD (IntegerField, ~2 billion limit) # BIGAUTOFIELD WILL BE REQUIRED FOR MORE! - - survey = models.PositiveSmallIntegerField() - - hpid = models.IntegerField(db_index=True) # up to over 200 million - - counts = models.IntegerField() # f4, up to ~44k integer: 2 byte too small - - exposure = models.FloatField() # f4, up to ~13300 float - + + survey = models.PositiveSmallIntegerField() + + hpid = models.IntegerField(db_index=True) # up to over 200 million + + counts = models.IntegerField() # f4, up to ~44k integer: 2 byte too small + + exposure = models.FloatField() # f4, up to ~13300 float + contaminated = models.BooleanField(default=False) - + def __str__(self): return f"Pixel {self.id} hpid {self.hpid} (Survey {self.survey.number})" - - class CatalogSource(models.Model): - - srcid = models.SmallIntegerField(primary_key=True) - - name = models.CharField(max_length=21) - - ra_deg = models.FloatField() - dec_deg = models.FloatField() - - pos_error = models.FloatField() - + srcid = models.SmallIntegerField(primary_key=True) + + name = models.CharField(max_length=21) + + ra_deg = models.FloatField() + + dec_deg = models.FloatField() + + pos_error = models.FloatField() + significance = models.FloatField() - - flux = models.FloatField() - - flux_error = models.FloatField() - + + flux = models.FloatField() + + flux_error = models.FloatField() + catalog_name = models.CharField(max_length=28) - - new_xray = models.BooleanField(default=False) - - source_type = models.CharField(max_length=13) \ No newline at end of file + + new_xray = models.BooleanField(default=False) + + source_type = models.CharField(max_length=13) diff --git a/serializers.py b/serializers.py index 3a68ff1..99430c8 100644 --- a/serializers.py +++ b/serializers.py @@ -1,8 +1,3 @@ # uplim/serializers.py from rest_framework import serializers from uplim.models import Pixel - -# class PixelSerializer(serializers.ModelSerializer): -# class Meta: -# model = Pixel -# fields = ('hpid', 'counts', 'exposure', 'contaminated') diff --git a/urls.py b/urls.py index ac44885..7502b60 100644 --- a/urls.py +++ b/urls.py @@ -1,10 +1,10 @@ # uplim/urls.py from django.urls import path -from .views import PixelAggregateView, UpperLimitView #, PixelDetailView +from .views import PixelAggregateView, UpperLimitView # , PixelDetailView urlpatterns = [ - #path('pixel//', PixelDetailView.as_view(), name='pixel-detail'), - path('pixel-aggregate/', PixelAggregateView.as_view(), name='pixel-aggregate'), - path('upper-limit/', UpperLimitView.as_view(), name='upper-limit'), + # path('pixel//', PixelDetailView.as_view(), name='pixel-detail'), + path("pixel-aggregate/", PixelAggregateView.as_view(), name="pixel-aggregate"), + path("upper-limit/", UpperLimitView.as_view(), name="upper-limit"), ] diff --git a/views.py b/views.py index 68b5eba..7d6dc3b 100644 --- a/views.py +++ b/views.py @@ -4,25 +4,24 @@ # search for pixels non-inclusively import healpy as hp import astropy.units as u -from astropy.coordinates import SkyCoord, Angle import numpy as np - import scipy.special as sp + +from astropy.coordinates import SkyCoord, Angle from astropy.stats import poisson_conf_interval from django.db.models import Sum +from django.shortcuts import get_object_or_404 from rest_framework.views import APIView from rest_framework.response import Response from rest_framework import status - -from django.shortcuts import get_object_or_404 - from uplim.models import Pixel, CatalogSource # SANITIZE RESPONSE DATA BEFORE JSON CONVERSION FOR DEBUGGING NANS # now NaNs are converted to 'null' beforehand # **************************************************************** + def sanitize(obj): if isinstance(obj, dict): return {k: sanitize(v) for k, v in obj.items()} @@ -52,17 +51,17 @@ def parse_survey_param(raw): return sorted(surveys) - # PIXEL VIEW (MOSTLY FOR TESTING) # add healpix indices into the output # ************************************************************** + class PixelAggregateView(APIView): def get(self, request): # GET PARAMETERS FROM THE QUERY # ************************************************************** - raw_pixel = request.query_params.get("pixel") + raw_pixel = request.query_params.get("pixel") raw_survey = request.query_params.get("survey") # 400 BADREQUEST @@ -70,7 +69,7 @@ class PixelAggregateView(APIView): if raw_pixel is None or raw_survey is None: return Response( {"detail": "Both `pixel` and `survey` parameters are required."}, - status=status.HTTP_400_BAD_REQUEST + status=status.HTTP_400_BAD_REQUEST, ) # FILTER THE INPUTS @@ -80,7 +79,7 @@ class PixelAggregateView(APIView): except ValueError: return Response( {"detail": "`pixel` must be an integer."}, - status=status.HTTP_400_BAD_REQUEST + status=status.HTTP_400_BAD_REQUEST, ) try: @@ -88,44 +87,33 @@ class PixelAggregateView(APIView): except ValueError: return Response( {"detail": "Malformed `survey`; use N, N,M or N-M formats."}, - status=status.HTTP_400_BAD_REQUEST + status=status.HTTP_400_BAD_REQUEST, ) # FILTER AND AGGREGATE # ************************************************************** - qs = Pixel.objects.filter( - hpid=hpid, - survey__in=survey_numbers - ) - + qs = Pixel.objects.filter(hpid=hpid, survey__in=survey_numbers) + if not qs.exists(): # no matching pixel(s) → 404 - get_object_or_404( - Pixel, - hpid=hpid, - survey__in=survey_numbers - ) + get_object_or_404(Pixel, hpid=hpid, survey__in=survey_numbers) aggregates = qs.aggregate( - #pixel_hpid=hpid, - #survey_number=survey_numbers, + # pixel_hpid=hpid, + # survey_number=survey_numbers, total_counts=Sum("counts"), - total_exposure=Sum("exposure") + total_exposure=Sum("exposure"), ) - plusdata = { - 'pixel_hpid' : hpid, - 'surveys' : survey_numbers - } - + plusdata = {"pixel_hpid": hpid, "surveys": survey_numbers} + result = {**aggregates, **plusdata} - + # RETURN THE SUMS # ************************************************************** return Response(result, status=status.HTTP_200_OK) - # UPPER LIMIT COMPUTATION VIEW # ************************************************************** @@ -134,317 +122,303 @@ class UpperLimitView(APIView): """ Calculate confidence bounds based on aperture photometry using classic and bayesian methods """ + def get(self, request): - + # GET PARAMETERS FROM THE REQUEST # ************************************************************** - + try: - ra = float(request.query_params.get('ra')) - dec = float(request.query_params.get('dec')) - confidence_level = float(request.query_params.get('cl')) + ra = float(request.query_params.get("ra")) + dec = float(request.query_params.get("dec")) + confidence_level = float(request.query_params.get("cl")) except (TypeError, ValueError): return Response( {"error": "Invalud parameters, provide RA, DEC, and CL"}, - status = status.HTTP_400_BAD_REQUEST + status=status.HTTP_400_BAD_REQUEST, ) # ── NEW: pull & parse survey selection ── - raw_survey = request.query_params.get('survey') + raw_survey = request.query_params.get("survey") if raw_survey is None: return Response( {"error": "Missing required `survey` parameter (e.g. ?survey=1,3-5)"}, - status=status.HTTP_400_BAD_REQUEST + status=status.HTTP_400_BAD_REQUEST, ) try: survey_numbers = parse_survey_param(raw_survey) except ValueError: return Response( - {"error": "Malformed `survey`; use formats like 1, 2-5, or 1,3-4"}, - status=status.HTTP_400_BAD_REQUEST - ) + {"error": "Malformed `survey`; use formats like 1, 2-5, or 1,3-4"}, + status=status.HTTP_400_BAD_REQUEST, + ) # hp = HEALPix( # nside = 4096, # order = 'ring', # frame = 'icrs' # ) - + # CREATE SKYCOORD, CONVERT TO GALACTIC BECAUSE THE HEALPIX # MAP ITSELF WAS MADE IN GALACTIC COORDINATES # ************************************************************** - - src_coord = SkyCoord( - ra, dec, unit = 'deg', frame = 'icrs' - ) + + src_coord = SkyCoord(ra, dec, unit="deg", frame="icrs") gal = src_coord.galactic - src_vec = hp.ang2vec(gal.l.deg, gal.b.deg, lonlat = True) + src_vec = hp.ang2vec(gal.l.deg, gal.b.deg, lonlat=True) # DEFINE APERTURE AND ANNULUS RADII # ************************************************************** - aperture_radius = 71 # radius of the aperture in arc seconds + aperture_radius = 71 # radius of the aperture in arc seconds # HPD ~48 arcseconds # 90% ~100 arcseconds - annulus_inner = 142 # 2 * aperture_radius - annulus_outer = 284 # 4 * aperture_radius - + annulus_inner = 142 # 2 * aperture_radius + annulus_outer = 284 # 4 * aperture_radius + # FETCH PIXEL DATA DEFINED VIA HP.QUERY_DISC (INCLUSIVE=FALSE) # ************************************************************** - + source_pixel_list = hp.query_disc( - nside = 4096, - vec = src_vec, - inclusive = False, - nest = False, - radius = (aperture_radius * u.arcsecond).to(u.radian).value + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(aperture_radius * u.arcsecond).to(u.radian).value, ) - + inner_pixel_list = hp.query_disc( - nside = 4096, - vec = src_vec, - inclusive = False, - nest = False, - radius = (annulus_inner * u.arcsecond).to(u.radian).value + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(annulus_inner * u.arcsecond).to(u.radian).value, ) - + outer_pixel_list = hp.query_disc( - nside = 4096, - vec = src_vec, - inclusive = False, - nest = False, - radius = (annulus_outer * u.arcsecond).to(u.radian).value + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(annulus_outer * u.arcsecond).to(u.radian).value, ) - - + annulus_pixel_list = [ item for item in outer_pixel_list if item not in inner_pixel_list ] - - + source_pixels = Pixel.objects.filter( - hpid__in = source_pixel_list, - survey__in = survey_numbers - ) - - annulus_pixels = Pixel.objects.filter( - hpid__in = annulus_pixel_list, - survey__in = survey_numbers - ) - - # check contamination - contamination = ( - source_pixels.filter(contaminated=True).exists() or - annulus_pixels.filter(contaminated=True).exists() + hpid__in=source_pixel_list, survey__in=survey_numbers + ) + + annulus_pixels = Pixel.objects.filter( + hpid__in=annulus_pixel_list, survey__in=survey_numbers + ) + + # check contamination + contamination = ( + source_pixels.filter(contaminated=True).exists() + or annulus_pixels.filter(contaminated=True).exists() ) - if not source_pixels.exists() and not annulus_pixels.exists(): return Response( {"detail": "No pixel data for the given survey selection."}, - status=status.HTTP_404_NOT_FOUND - ) + status=status.HTTP_404_NOT_FOUND, + ) - - # COMPUTE COUNTS, BACKGROUND ESTIMATE, EXPOSURE # ************************************************************** - + N = sum(obj.counts for obj in source_pixels) - + Nnpix = len(source_pixels) - + Bcounts = sum(obj.counts for obj in annulus_pixels) - + Bnpix = len(annulus_pixels) - + B = Bcounts / Bnpix * Nnpix - + tsum = sum(obj.exposure for obj in source_pixels) - + t = tsum / Nnpix - + # CONSTANTS # ************************************************************** - - #EEF = .9 # eclosed energy fraction, .5 for hpd, .9 for w90 - #ECF = 4e-11 # energy conversion factor - - EEF = .80091 # use values from the paper + + # EEF = .9 # eclosed energy fraction, .5 for hpd, .9 for w90 + # ECF = 4e-11 # energy conversion factor + + EEF = 0.80091 # use values from the paper ECF = 3.3423184e-11 - # BAYESIAN IMPLEMENTATION VIA POISSON_CONF_INTERVAL # ************************************************************** - + low, high = poisson_conf_interval( - n = N, - background = B, - interval = 'kraft-burrows-nousek', - confidence_level=confidence_level + n=N, + background=B, + interval="kraft-burrows-nousek", + confidence_level=confidence_level, ) - + bayesian_count_ul = high bayesian_count_ll = low - - bayesian_rate_ul = bayesian_count_ul / t / EEF # count rate limits + + bayesian_rate_ul = bayesian_count_ul / t / EEF # count rate limits bayesian_rate_ll = bayesian_count_ll / t / EEF - - bayesian_flux_ul = bayesian_rate_ul * ECF # flux limits + + bayesian_flux_ul = bayesian_rate_ul * ECF # flux limits bayesian_flux_ll = bayesian_rate_ll * ECF - + # CLASSICAL IMPLEMENTATION VIA GAMMAINCCINV # **************************************************************** - - classic_count_ul = sp.gammainccinv(N+1, 1 - confidence_level) - B + + classic_count_ul = sp.gammainccinv(N + 1, 1 - confidence_level) - B classic_count_ll = sp.gammainccinv(N, confidence_level) - B - + if not np.isfinite(classic_count_ll) or classic_count_ll < 0: classic_count_ll = 0.0 - - - classic_rate_ul = classic_count_ul / t / EEF # count rate limits - classic_rate_ll = classic_count_ll / t / EEF - - classic_flux_ul = classic_rate_ul * ECF # flux limits + + classic_rate_ul = classic_count_ul / t / EEF # count rate limits + classic_rate_ll = classic_count_ll / t / EEF + + classic_flux_ul = classic_rate_ul * ECF # flux limits classic_flux_ll = classic_rate_ll * ECF - + # FLUX ESTIMATION # **************************************************************** - - S = N - B # counts as simply counts within aperture - # with the background estimate subtracted - - - CR = S / t / EEF # source count rate - - BR = B / t # background rate within aperture - - FL = CR * ECF # conversion to flux - - Flux = max(FL, 0) # flux cannot be lower than zero - + + S = N - B # counts as simply counts within aperture + # with the background estimate subtracted + + CR = S / t / EEF # source count rate + + BR = B / t # background rate within aperture + + FL = CR * ECF # conversion to flux + + Flux = max(FL, 0) # flux cannot be lower than zero + # NEARBY SOURCES CHECK # **************************************************************** - + radius_as = 120 radius_deg = radius_as / 3600 - + dec_min = max(dec - radius_deg, -90) - dec_max = min(dec + radius_deg, 90) - + dec_max = min(dec + radius_deg, 90) + # cheap belt query belt_sources = CatalogSource.objects.filter( - dec_deg__gte = dec_min, - dec_deg__lte = dec_max + dec_deg__gte=dec_min, dec_deg__lte=dec_max ) - - center_coord = SkyCoord(ra, dec, unit='deg') - + + center_coord = SkyCoord(ra, dec, unit="deg") + nearby_sources = [] - - #refine belt to circular region using astropy separation + + # refine belt to circular region using astropy separation for catsrc in belt_sources: - catsrc_coord = SkyCoord(catsrc.ra_deg, catsrc.dec_deg, unit='deg') + catsrc_coord = SkyCoord(catsrc.ra_deg, catsrc.dec_deg, unit="deg") if center_coord.separation(catsrc_coord).deg <= radius_deg: nearby_sources.append( { - 'srcid' : catsrc.srcid, - 'name' : catsrc.name, - 'ra_deg' : catsrc.ra_deg, - 'dec_deg' : catsrc.dec_deg, - 'pos_error' : catsrc.pos_error, - 'significance' : catsrc.significance, - 'flux' : catsrc.flux, - 'flux_error' : catsrc.flux_error, - 'catalog_name' : catsrc.catalog_name, - 'new_xray' : catsrc.new_xray, - 'source_type' : catsrc.source_type + "srcid": catsrc.srcid, + "name": catsrc.name, + "ra_deg": catsrc.ra_deg, + "dec_deg": catsrc.dec_deg, + "pos_error": catsrc.pos_error, + "significance": catsrc.significance, + "flux": catsrc.flux, + "flux_error": catsrc.flux_error, + "catalog_name": catsrc.catalog_name, + "new_xray": catsrc.new_xray, + "source_type": catsrc.source_type, } ) - # SQUARE REGION IMAGE SERVING # **************************************************************** - - # get hpids within a circle of radius sqrt(2) * outer annulus radius + + # get hpids within a circle of radius sqrt(2) * outer annulus radius map_radius = annulus_outer * np.sqrt(2) - + map_pixel_list = hp.query_disc( - nside = 4096, - vec = src_vec, - inclusive = False, - nest = False, - radius = (map_radius * u.arcsecond).to(u.radian).value + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(map_radius * u.arcsecond).to(u.radian).value, ) - + # fetch those pixels for the requested surveys # summing counts and sorting by hpid - + map_pixels_qs = ( - Pixel.objects - .filter(hpid__in = map_pixel_list, survey__in = survey_numbers) - .values('hpid') - .annotate(counts=Sum('counts')) - .order_by('hpid') + Pixel.objects.filter(hpid__in=map_pixel_list, survey__in=survey_numbers) + .values("hpid") + .annotate(counts=Sum("counts")) + .order_by("hpid") ) - + # turn the queryset to a list map_pixels_list = list(map_pixels_qs) - + # get lists of healpix indices and count values - map_healpix_list = [d['hpid'] for d in map_pixels_list] - map_counts_list = [d['counts'] for d in map_pixels_list] - + map_healpix_list = [d["hpid"] for d in map_pixels_list] + map_counts_list = [d["counts"] for d in map_pixels_list] + # set map nside map_nside = 4096 - + # set map order - map_order = 'ring' - + map_order = "ring" + # assemble the result dict map_dict = { - 'healpix' : map_healpix_list, - 'counts' : map_counts_list, - 'nside' : map_nside, - 'order' : map_order, - 'radius_as' : map_radius + "healpix": map_healpix_list, + "counts": map_counts_list, + "nside": map_nside, + "order": map_order, + "radius_as": map_radius, } - + # RESULT JSON # **************************************************************** - + result = { # frequentist limits - 'ClassicUpperLimit' : classic_count_ul, - 'ClassicLowerLimit' : classic_count_ll, - 'ClassicCountRateUpperLimit' : classic_rate_ul, - 'ClassicCountRateLowerLimit' : classic_rate_ll, - 'ClassicFluxUpperLimit' : classic_flux_ul, - 'ClassicFluxLowerLimit' : classic_flux_ll, + "ClassicUpperLimit": classic_count_ul, + "ClassicLowerLimit": classic_count_ll, + "ClassicCountRateUpperLimit": classic_rate_ul, + "ClassicCountRateLowerLimit": classic_rate_ll, + "ClassicFluxUpperLimit": classic_flux_ul, + "ClassicFluxLowerLimit": classic_flux_ll, # bayesian limits - 'BayesianUpperLimit' : bayesian_count_ul, - 'BayesianLowerLimit' : bayesian_count_ll, - 'BayesianCountRateUpperLimit' : bayesian_rate_ul, - 'BayesianCountRateLowerLimit' : bayesian_rate_ll, - 'BayesianFluxUpperLimit' : bayesian_flux_ul, - 'BayesianFluxLowerLimit' : bayesian_flux_ll, + "BayesianUpperLimit": bayesian_count_ul, + "BayesianLowerLimit": bayesian_count_ll, + "BayesianCountRateUpperLimit": bayesian_rate_ul, + "BayesianCountRateLowerLimit": bayesian_rate_ll, + "BayesianFluxUpperLimit": bayesian_flux_ul, + "BayesianFluxLowerLimit": bayesian_flux_ll, # flux 'center value' estimate - 'FluxEstimate' : Flux, + "FluxEstimate": Flux, # raw data - 'ApertureCounts' : N, - 'ApertureBackgroundCounts' : B, - 'SourceCounts' : S, - 'Exposure' : t, + "ApertureCounts": N, + "ApertureBackgroundCounts": B, + "SourceCounts": S, + "Exposure": t, # count rates - 'SourceRate' : CR, - 'BackgroundRate' : BR, + "SourceRate": CR, + "BackgroundRate": BR, # contamination - 'Contamination' : contamination, - 'NearbySources' : nearby_sources, + "Contamination": contamination, + "NearbySources": nearby_sources, # count map for the frontend image - 'CountMap' : map_dict + "CountMap": map_dict, } - - clean = sanitize(result) # calling sanitize() to convert NaN to null - - return Response(clean, status=status.HTTP_200_OK) \ No newline at end of file + + clean = sanitize(result) # calling sanitize() to convert NaN to null + + return Response(clean, status=status.HTTP_200_OK)