initial commit

admin.py

@@ -0,0 +1,3 @@
+from django.contrib import admin
+
+# Register your models here.

apps.py

@@ -0,0 +1,6 @@
+from django.apps import AppConfig
+
+
+class AxcUlConfig(AppConfig):
+    default_auto_field = 'django.db.models.BigAutoField'
+    name = 'axc_ul'

management/commands/load_survey.py

@@ -0,0 +1,110 @@
+import numpy as np
+from astropy.io import fits
+
+from django.core.management.base import BaseCommand
+from django.db import transaction
+from axc_ul.models import Pixel
+
+from itertools import islice
+
+BATCH_SIZE = 1000000
+
+def batch(iterable, size):
+    """
+    Generator that yields successive chunks of size 'size' from 'iterable'.
+    """
+    iterable = iter(iterable)
+    while True:
+        chunk = list(islice(iterable, size))
+        if not chunk:
+            break
+        yield chunk
+
+
+class Command(BaseCommand):
+    help = "Process FITS files and store the data in the database"
+
+    def add_arguments(self, parser):
+        parser.add_argument(
+            '--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'
+        )
+        # parser.add_argument(
+        #     '--survey_number', 
+        #     type=int, 
+        #     required=True,
+        #     help='Integer ID of the survey being read'
+        # )
+
+    def handle(self, *args, **options):
+        
+        counts_file = options['counts']
+        exposure_file = options['exposure']
+        # survey_number = options['survey_number']
+
+        self.stdout.write(f"Counts file:   {counts_file}")
+        self.stdout.write(f"Exposure file: {exposure_file}")
+
+        with fits.open(counts_file) as hdul:
+
+            column_name = "T"
+            counts_map = hdul[1].data[column_name]
+
+        counts_data = counts_map.ravel()
+
+
+        with fits.open(exposure_file) as hdul:
+
+            column_name = "T"
+            exposure_map = hdul[1].data[column_name]
+
+        exposure_data = exposure_map.ravel()
+
+
+
+        self.stdout.write(f"Counts Data Shape:   {counts_data.shape}")
+        self.stdout.write(f"Exposure Data Shape: {exposure_data.shape}")
+
+        assert counts_data.shape == exposure_data.shape, "Counts and exposure maps must have the same shape"
+
+        #rate_data = np.divide(counts_data, exposure_data)
+        
+        # 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}")
+        
+
+        # Create a generator that yields Pixel objects one by one.
+        pixel_generator = (
+            Pixel(
+                hpid=i,
+                counts=int(count),
+                exposure=float(exposure),
+                #rate=float(rate),
+                #survey=survey
+            )
+            for i, (count, exposure) in enumerate(zip(counts_data, exposure_data))
+        )
+
+        total_inserted = 0
+        # Process the generator in batches.
+        for pixel_batch in batch(pixel_generator, BATCH_SIZE):
+            with transaction.atomic():
+                Pixel.objects.bulk_create(pixel_batch)
+            total_inserted += len(pixel_batch)
+            self.stdout.write(f"Inserted {total_inserted} pixels")
+
+        self.stdout.write(f"Inserted a total of {total_inserted} pixels.")

migrations/0001_initial.py

@@ -0,0 +1,35 @@
+# Generated by Django 5.1.6 on 2025-02-27 12:55
+
+import django.db.models.deletion
+import uuid
+from django.db import migrations, models
+
+
+class Migration(migrations.Migration):
+
+    initial = True
+
+    dependencies = [
+    ]
+
+    operations = [
+        migrations.CreateModel(
+            name='Survey',
+            fields=[
+                ('id', models.BigAutoField(auto_created=True, primary_key=True, serialize=False, verbose_name='ID')),
+                ('number', models.IntegerField()),
+            ],
+        ),
+        migrations.CreateModel(
+            name='Pixel',
+            fields=[
+                ('id', models.BigAutoField(auto_created=True, primary_key=True, serialize=False, verbose_name='ID')),
+                ('uuid', models.UUIDField(default=uuid.uuid4, editable=False, unique=True)),
+                ('hpid', models.IntegerField()),
+                ('counts', models.IntegerField()),
+                ('exposure', models.FloatField()),
+                ('rate', models.FloatField()),
+                ('survey', models.ForeignKey(on_delete=django.db.models.deletion.CASCADE, related_name='pixels', to='axc_ul.survey')),
+            ],
+        ),
+    ]

migrations/0002_remove_pixel_uuid_alter_pixel_id.py

@@ -0,0 +1,22 @@
+# Generated by Django 5.1.6 on 2025-02-27 19:35
+
+from django.db import migrations, models
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('axc_ul', '0001_initial'),
+    ]
+
+    operations = [
+        migrations.RemoveField(
+            model_name='pixel',
+            name='uuid',
+        ),
+        migrations.AlterField(
+            model_name='pixel',
+            name='id',
+            field=models.AutoField(primary_key=True, serialize=False),
+        ),
+    ]

migrations/0003_alter_pixel_hpid.py

@@ -0,0 +1,18 @@
+# Generated by Django 5.1.6 on 2025-03-02 15:09
+
+from django.db import migrations, models
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('axc_ul', '0002_remove_pixel_uuid_alter_pixel_id'),
+    ]
+
+    operations = [
+        migrations.AlterField(
+            model_name='pixel',
+            name='hpid',
+            field=models.IntegerField(db_index=True),
+        ),
+    ]

migrations/0004_survey_nside.py

@@ -0,0 +1,18 @@
+# Generated by Django 5.1.6 on 2025-03-02 16:37
+
+from django.db import migrations, models
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('axc_ul', '0003_alter_pixel_hpid'),
+    ]
+
+    operations = [
+        migrations.AddField(
+            model_name='survey',
+            name='nside',
+            field=models.IntegerField(default=4096),
+        ),
+    ]

migrations/0005_remove_pixel_survey_remove_pixel_rate_delete_survey.py

@@ -0,0 +1,24 @@
+# Generated by Django 5.1.6 on 2025-03-06 08:14
+
+from django.db import migrations
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('axc_ul', '0004_survey_nside'),
+    ]
+
+    operations = [
+        migrations.RemoveField(
+            model_name='pixel',
+            name='survey',
+        ),
+        migrations.RemoveField(
+            model_name='pixel',
+            name='rate',
+        ),
+        migrations.DeleteModel(
+            name='Survey',
+        ),
+    ]

migrations/0006_pixel_contaminated.py

@@ -0,0 +1,18 @@
+# Generated by Django 5.1.6 on 2025-03-13 13:27
+
+from django.db import migrations, models
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('axc_ul', '0005_remove_pixel_survey_remove_pixel_rate_delete_survey'),
+    ]
+
+    operations = [
+        migrations.AddField(
+            model_name='pixel',
+            name='contaminated',
+            field=models.BooleanField(default=False),
+        ),
+    ]

models.py

@@ -0,0 +1,14 @@
+from django.db import models
+
+class Pixel(models.Model):
+    
+    id = models.AutoField(primary_key=True)
+    
+    hpid = models.IntegerField(db_index=True)
+    
+    counts = models.IntegerField()
+    
+    exposure = models.FloatField()
+    
+    contaminated = models.BooleanField(default=False)
+    

serializers.py

@@ -0,0 +1,8 @@
+# serializers.py
+from rest_framework import serializers
+from axc_ul.models import Pixel
+
+class PixelSerializer(serializers.ModelSerializer):
+    class Meta:
+        model = Pixel
+        fields = ('hpid', 'counts', 'exposure', 'contaminated')

tests.py

@@ -0,0 +1,3 @@
+from django.test import TestCase
+
+# Create your tests here.

urls.py

@@ -0,0 +1,8 @@
+# urls.py
+from django.urls import path
+from .views import PixelDetailView, UpperLimitView
+
+urlpatterns = [
+    path('pixel/<int:hpid>/', PixelDetailView.as_view(), name='pixel-detail'),
+    path('upper-limit/', UpperLimitView.as_view(), name='upper-limit')
+]

views.py

@@ -0,0 +1,157 @@
+# axc_ul/views.py
+
+# from astropy_healpix import HEALPix does not have an option to
+#                                     search for pixels non-inclusively
+import healpy as hp
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+import scipy.special as sp
+
+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 axc_ul.models import Pixel
+
+from .serializers import PixelSerializer
+
+class PixelDetailView(APIView):
+    """
+    API endpoint that returns the pixel data (counts, exposure, rate)
+    for a given survey number and hpid.
+    """
+    def get(self, request, hpid):
+        # Get the survey using the survey_number field.
+        # survey = get_object_or_404(Survey, number=survey_number)
+        
+        # Get the pixel corresponding to the survey and hpid.
+        pixel = get_object_or_404(Pixel, hpid=hpid)
+        
+        # Serialize the pixel data to JSON.
+        serializer = PixelSerializer(pixel)
+        return Response(serializer.data, status=status.HTTP_200_OK)
+
+
+class UpperLimitView(APIView):
+    """
+    Calculate and return upper limits (UL, CRUL, FLUL) based on provided RA, DEC, and CL parameters.
+
+    Args:
+        request: HTTP GET request containing 'ra', 'dec', and 'cl' query parameters.
+
+    Returns:
+        Response object with JSON data containing calculated UL, CRUL, and FLUL values.
+
+    Raises:
+        ValueError if provided parameters are invalid.
+    """
+    def get(self, request):
+        
+        try:
+            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
+            )
+        
+        # hp = HEALPix(
+        #     nside = 4096,
+        #     order = 'ring',
+        #     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)
+
+        aperture_radius = 50 # radius of the aperture in arc seconds
+        # HPD ~48 arcseconds
+        # 90% ~100 arcseconds
+        annulus_inner = 100 # 2 * aperture_radius
+        annulus_outer = 200 # 4 * aperture_radius
+        
+        source_pixel_list = hp.query_disc(
+            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
+        )
+        
+        outer_pixel_list = hp.query_disc(
+            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) 
+        
+        annulus_pixels = Pixel.objects.filter(hpid__in = annulus_pixel_list)
+        
+        
+        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
+        
+        
+        EEF = .9      # eclosed energy fraction, .5 for hpd, .9 for w90
+        ECF = 4e-11   # energy conversion factor
+        
+        UL = (
+            sp.gammaincinv(
+                N + 1, 
+                confidence_level * sp.gammaincc(N + 1, B) + sp.gammainc(N + 1, B)
+            ) - B
+        )                   # count upper limit
+                            
+         
+        CRUL = UL / t / EEF # count rate upper limit
+        
+        FLUL = CRUL * ECF   # flux upper limit
+        
+        return Response({
+            'UpperLimit' : UL,
+            'CountRateUpperLimit' : CRUL,
+            'FluxUpperLimit' : FLUL,
+            'N' : N,
+            'Nnpix' : Nnpix,
+            'Bcounts' : Bcounts,
+            'Bnpix' : Bnpix,
+            'B' : B,
+            'tsum' : tsum,
+            't' : t
+        }, status=status.HTTP_200_OK)