Merge pull request #145 from Yong2Sheng/fast_ts_map_galactic

Memory leak, downloading cells and SMEX GRB
cositools · Feb 28, 2024 · a109ad8 · a109ad8
2 parents bfaa42f + c720d41
commit a109ad8
Show file tree

Hide file tree

Showing 16 changed files with 533 additions and 286 deletions.
diff --git a/cosipy/ts_map/__init__.py b/cosipy/ts_map/__init__.py
@@ -1,2 +1,2 @@
 from .TSMap import TSMap
-from .fast_ts_fit import FastTSMap
+from .fast_ts_fit import FastTSMap
diff --git a/cosipy/ts_map/fast_ts_fit.py b/cosipy/ts_map/fast_ts_fit.py
@@ -14,6 +14,9 @@
 from cosipy.response import FullDetectorResponse
 import time
 import scipy.stats
+import os
+import psutil
+import gc
 
 class FastTSMap():
 
@@ -134,6 +137,10 @@ def get_cds_array(hist, energy_channel):
         hist_cds_sliced = FastTSMap.slice_energy_channel(hist, energy_channel[0], energy_channel[1])   
         hist_cds = hist_cds_sliced.project(["PsiChi", "Phi"])
         cds_array = np.array(hist_cds.to_dense()[:]).flatten()  # here [:] is equivalent to [:, :]
+        del hist
+        del hist_cds_sliced
+        del hist_cds
+        gc.collect()
 
         return cds_array
 
@@ -216,23 +223,37 @@ def get_ei_cds_array(hypothesis_coord, energy_channel, response_path, spectrum,
         """
 
         # check inputs, will complete later
-
-        # the local and galactic frame works very differently, so we need to compuate the point source response (psr) accordingly 
+
+        # the local and galactic frame works very differently, so we need to compuate the point source response (psr) accordingly
+        #time_cds_start = time.time()
         if cds_frame == "local":
 
             if orientation == None:
                 raise TypeError("The when the data are binned in local frame, orientation must be provided to compute the expected counts.")
-
+
+            #time_coord_convert_start = time.time()
             # convert the hypothesis coord to the local frame (Spacecraft frame)
             hypothesis_in_sc_frame = orientation.get_target_in_sc_frame(target_name = "Hypothesis", 
                                                                         target_coord = hypothesis_coord, 
                                                                         quiet = True)
+            #time_coord_convert_end = time.time()
+            #time_coord_convert_used = time_coord_convert_end - time_coord_convert_start
+            #print(f"The time used for coordinate conversion is {time_coord_convert_used}s.")
+
+            #time_dwell_start = time.time()
             # get the dwell time map: the map of the time spent on each pixel in the local frame
             dwell_time_map = orientation.get_dwell_map(response = response_path)
+            #time_dwell_end = time.time()
+            #time_dwell_used = time_dwell_end - time_dwell_start
+            #print(f"The time used for dwell time map is {time_dwell_used}s.")
 
+            #time_psr_start = time.time()
             # convolve the response with the dwell_time_map to get the point source response
             with FullDetectorResponse.open(response_path) as response:
                 psr = response.get_point_source_response(dwell_time_map)
+            #time_psr_end = time.time()
+            #time_psr_used = time_psr_end - time_psr_start
+            #print(f"The time used for psr is {time_psr_used}s.")
 
         elif cds_frame == "galactic":
 
@@ -243,9 +264,17 @@ def get_ei_cds_array(hypothesis_coord, energy_channel, response_path, spectrum,
 
         # convolve the point source reponse with the spectrum to get the expected counts
         expectation = psr.get_expectation(spectrum)
+        del psr
+        gc.collect()
 
         # slice energy channals and project it to CDS
         ei_cds_array = FastTSMap.get_cds_array(expectation, energy_channel)
+        del expectation
+        gc.collect()
+
+        #time_cds_end = time.time()
+        #time_cds_used = time_cds_end - time_cds_start
+        #print(f"The time used for cds is {time_cds_used}s.")
 
         return ei_cds_array
 
@@ -439,7 +468,7 @@ def _plot_ts(result_array, skycoord = None, containment = None):
 
         return
 
-    def plot_ts(self, skycoord = None, containment = 0.9):
+    def plot_ts(self, skycoord = None, containment = None):
 
         """
         Plot the containment region of the TS map.
@@ -481,6 +510,14 @@ def get_chi_critical_value(containment = 0.90):
         """
 
         return scipy.stats.chi2.ppf(containment, df=2)
+
+    @staticmethod
+    def show_memory_info(hint):
+        pid = os.getpid()
+        p = psutil.Process(pid)
 
+        info = p.memory_full_info()
+        memory = info.uss / 1024. / 1024
+        print('{} memory used: {} MB'.format(hint, memory))
 
 
diff --git a/docs/tutorials/CygX1_source_path_in_SC_frame.npy b/docs/tutorials/CygX1_source_path_in_SC_frame.npy
diff --git a/docs/tutorials/Intro.ipynb b/docs/tutorials/Intro.ipynb
diff --git a/docs/tutorials/data/background.yaml b/docs/tutorials/data/background.yaml
diff --git a/docs/tutorials/data/grb.yaml b/docs/tutorials/data/grb.yaml
diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst
@@ -14,13 +14,14 @@ List of tutorials and contents, as a link to the corresponding Python notebook i
   - Combining files.
   - Inspecting and plotting the data
 
-2. Spacecraft orientation and location `(ipynb) <https://github.com/cositools/cosipy/tree/main/docs/tutorials/Point_source_resonse.ipynb>`_
+2. Spacecraft orientation and location `(ipynb) <https://github.com/cositools/cosipy/tree/main/docs/tutorials/response/Point_source_response.ipynb>`_
 
   - SC file format and manipulation it —e.g. get a time range, rebin it.  
   - The dwell time map and how to obtain it
+  - Generate point source response and export to the format that can be read by XSPEC
   - The scatt map and how to obtain it
 
-3. Detector response and signal expectation `(ipynb) <https://github.com/cositools/cosipy/tree/main/docs/tutorials/DetectorResponse.ipynb>`_
+3. Detector response and signal expectation `(ipynb) <https://github.com/cositools/cosipy/tree/main/docs/tutorials/response/DetectorResponse.ipynb>`_
 
   - Explanation of the detector response format and meaning
   - Visualizing the response
@@ -55,7 +56,7 @@ List of tutorials and contents, as a link to the corresponding Python notebook i
   - Scatt binning and its advantages/disadvantages
   - Fitting the 511 diffuse emission.
 
-9. TODO: Source injector
+9. TODO: Source injector `(ipynb) <https://github.com/cositools/cosipy/tree/main/docs/tutorials/source_injector/GRB_source_injector.ipynb`
   - Nice to have: allow theorist to test the sensitivity of their models
 
 .. warning::

diff --git a/docs/tutorials/CygX1_DwellMap.fits → docs/tutorials/response/CygX1_DwellMap.fits b/docs/tutorials/CygX1_DwellMap.fits → docs/tutorials/response/CygX1_DwellMap.fits
diff --git a/docs/tutorials/DetectorResponse.ipynb → ...tutorials/response/DetectorResponse.ipynb b/docs/tutorials/DetectorResponse.ipynb → ...tutorials/response/DetectorResponse.ipynb
@@ -1271,9 +1271,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python [conda env:cosi_nomegalib]",
+   "display_name": "cosipy_spacecraftfile_new",
    "language": "python",
-   "name": "conda-env-cosi_nomegalib-py"
+   "name": "cosipy_spacecraftfile_new"
   },
   "language_info": {
    "codemirror_mode": {
@@ -1285,7 +1285,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.6"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,

diff --git a/docs/tutorials/Point_source_resonse.ipynb → ...ials/response/Point_source_response.ipynb b/docs/tutorials/Point_source_resonse.ipynb → ...ials/response/Point_source_response.ipynb
diff --git a/docs/tutorials/GRB_source_injector.ipynb → ...source_injector/GRB_source_injector.ipynb b/docs/tutorials/GRB_source_injector.ipynb → ...source_injector/GRB_source_injector.ipynb
diff --git a/docs/tutorials/ts_map/GRB_Orientation.ori b/docs/tutorials/ts_map/GRB_Orientation.ori