threeML / hawc_hal

HAWC Accelerated Likelihood - python-only framework for HAWC data analysis
BSD 3-Clause "New" or "Revised" License
11 stars 21 forks source link

Introducing multithreading for log-likelihood computation #93

Closed torresramiro350 closed 2 months ago

torresramiro350 commented 6 months ago
torresramiro350 commented 6 months ago

For comparison in speed, here's the execution time with and without multithreading. In this example, I fit a two extended source Gaussian_on_sphere with Log_parabola spectral shape. I fit the extension and spectral parameters for the Crab.

Without multithreading;

  _     ._   __/__   _ _  _  _ _/_   Recorded: 04:38:15  Samples:  930
 /_//_/// /_\ / //_// / //_'/ //     Duration: 276.934   CPU time: 260.048
/   _/                      v4.6.2

Program: crab_fitting.py -n 5

276.714 MainThread  <thread>:139920236996416
├─ 261.435 <module>  crab_fitting.py:1
│  └─ 261.435 main  crab_fitting.py:30
│     ├─ 198.165 JointLikelihood.fit  threeML/classicMLE/joint_likelihood.py:217
│     │     [8 frames hidden]  threeML
│     │        171.986 JointLikelihood.minus_log_like_profile  threeML/classicMLE/joint_likelihood.py:936
│     │        └─ 171.986 HAL.inner_fit  hawc_hal/HAL.py:1298
│     │           └─ 171.986 HAL.get_log_like  hawc_hal/HAL.py:858
│     │              └─ 171.236 HAL._get_expectation  hawc_hal/HAL.py:1001
│     │                 ├─ 135.947 ConvolvedExtendedSource2D.get_source_map  hawc_hal/convolved_source/convolved_extended_source.py:143
│     │                 │  ├─ 91.581 [self]  hawc_hal/convolved_source/convolved_extended_source.py
│     │                 │  ├─ 24.082 old_div  past/utils/__init__.py:85
│     │                 │  └─ 16.768 sum  <__array_function__ internals>:177
│     │                 │        [4 frames hidden]  <__array_function__ internals>, numpy...
│     │                 ├─ 25.524 PSFConvolutor.extended_source_image  hawc_hal/psf_fast/psf_convolutor.py:67
│     │                 │  ├─ 10.260 rfftn  <__array_function__ internals>:177
│     │                 │  │     [10 frames hidden]  <__array_function__ internals>, numpy...
│     │                 │  └─ 10.013 irfftn  <__array_function__ internals>:177
│     │                 │        [10 frames hidden]  <__array_function__ internals>, numpy...
│     │                 └─ 8.014 FlatSkyToHealpixTransform.__call__  hawc_hal/healpix_handling/flat_sky_to_healpix.py:101
│     │                    └─ 8.014 FastBilinearInterpolation.__call__  hawc_hal/interpolation/fast_bilinar_interpolation.py:88
│     │        24.580 JointLikelihood.minus_log_like_profile  threeML/classicMLE/joint_likelihood.py:936
│     │        └─ 24.330 HAL.inner_fit  hawc_hal/HAL.py:1298
│     │           └─ 24.330 HAL.get_log_like  hawc_hal/HAL.py:858
│     │              └─ 24.330 HAL._get_expectation  hawc_hal/HAL.py:1001
│     │                 ├─ 19.324 ConvolvedExtendedSource2D.get_source_map  hawc_hal/convolved_source/convolved_extended_source.py:143
│     │                 │  ├─ 11.545 [self]  hawc_hal/convolved_source/convolved_extended_source.py
│     │                 │  ├─ 4.020 old_div  past/utils/__init__.py:85
│     │                 │  └─ 3.256 sum  <__array_function__ internals>:177
│     │                 │        [4 frames hidden]  <__array_function__ internals>, numpy...
│     │                 └─ 3.754 PSFConvolutor.extended_source_image  hawc_hal/psf_fast/psf_convolutor.py:67
│     ├─ 32.967 HAL.__init__  hawc_hal/HAL.py:58
│     │  ├─ 19.060 map_tree_factory  hawc_hal/maptree/map_tree.py:21
│     │  │  └─ 19.060 MapTree.from_root_file  hawc_hal/maptree/map_tree.py:54
│     │  │     └─ 19.060 from_root_file  hawc_hal/maptree/from_root_file.py:165
│     │  │        └─ 19.060 Pool.starmap  multiprocessing/pool.py:369
│     │  │              [6 frames hidden]  multiprocessing, threading, <built-in>
│     │  └─ 13.155 hawc_response_factory  hawc_hal/response/response.py:31
│     │     └─ 13.155 HAWCResponse.from_root_file  hawc_hal/response/response.py:381
│     │        ├─ 9.000 ResponseBin.from_ttree  hawc_hal/response/response_bin.py:172
│     │        │  └─ 8.750 PSFWrapper.psf_eval  hawc_hal/psf_fast/psf_wrapper.py:222
│     │        │     └─ 5.250 PSFWrapper.__init__  hawc_hal/psf_fast/psf_wrapper.py:50
│     │        │        └─ 4.750 PSFWrapper._prepare_brightness_interpolation_points  hawc_hal/psf_fast/psf_wrapper.py:94
│     │        │           └─ 4.750 <listcomp>  hawc_hal/psf_fast/psf_wrapper.py:100
│     │        │              └─ 4.000 PSFWrapper.integral  hawc_hal/psf_fast/psf_wrapper.py:252
│     │        │                 └─ 3.750 InterpolatedUnivariateSpline.integral  scipy/interpolate/_fitpack2.py:423
│     │        │                       [2 frames hidden]  scipy
│     │        └─ 3.905 Pool.starmap  multiprocessing/pool.py:369
│     │              [6 frames hidden]  multiprocessing, threading, <built-in>
│     ├─ 13.615 JointLikelihood.compute_TS  threeML/classicMLE/joint_likelihood.py:1289
│     │     [2 frames hidden]  threeML
│     │        13.615 JointLikelihood._assign_model_to_data  threeML/classicMLE/joint_likelihood.py:118
│     │        └─ 13.615 HAL.set_model  hawc_hal/HAL.py:352
│     │           └─ 13.615 HAL._setup_psf_convolutors  hawc_hal/HAL.py:212
│     │              └─ 13.615 PSFConvolutor.__init__  hawc_hal/psf_fast/psf_convolutor.py:18
│     │                 └─ 13.615 PSFInterpolator.point_source_image  hawc_hal/psf_fast/psf_interpolator.py:123
│     │                    └─ 13.615 PSFInterpolator._get_point_source_image_aitoff  hawc_hal/psf_fast/psf_interpolator.py:36
│     │                       └─ 13.114 reproject_exact  reproject/spherical_intersect/high_level.py:11
│     │                             [4 frames hidden]  reproject
│     └─ 13.599 JointLikelihood.__init__  threeML/classicMLE/joint_likelihood.py:52
│           [2 frames hidden]  threeML
│              13.599 JointLikelihood._assign_model_to_data  threeML/classicMLE/joint_likelihood.py:118
│              └─ 13.599 HAL.set_model  hawc_hal/HAL.py:352
│                 └─ 13.599 HAL._setup_psf_convolutors  hawc_hal/HAL.py:212
│                    └─ 13.599 PSFConvolutor.__init__  hawc_hal/psf_fast/psf_convolutor.py:18
│                       └─ 13.599 PSFInterpolator.point_source_image  hawc_hal/psf_fast/psf_interpolator.py:123
│                          └─ 13.599 PSFInterpolator._get_point_source_image_aitoff  hawc_hal/psf_fast/psf_interpolator.py:36
│                             └─ 13.097 reproject_exact  reproject/spherical_intersect/high_level.py:11
│                                   [4 frames hidden]  reproject
└─ 13.279 log_likelihood  hawc_hal/log_likelihood.py:7

With multithreading:

  _     ._   __/__   _ _  _  _ _/_   Recorded: 04:32:51  Samples:  504
 /_//_/// /_\ / //_// / //_'/ //     Duration: 184.733   CPU time: 315.209
/   _/                      v4.6.2

Program: crab_fitting.py -n 5

184.517 <module>  crab_fitting.py:1
└─ 184.517 main  crab_fitting.py:30
   ├─ 119.175 JointLikelihood.fit  threeML/classicMLE/joint_likelihood.py:217
   │     [8 frames hidden]  threeML
   │        102.212 JointLikelihood.minus_log_like_profile  threeML/classicMLE/joint_likelihood.py:936
   │        └─ 101.962 HAL.inner_fit  hawc_hal/HAL.py:1609
   │           └─ 101.962 HAL.get_log_like  hawc_hal/HAL.py:1019
   │              └─ 101.962 HAL._multithreading_log_like  hawc_hal/HAL.py:974
   │                 └─ 101.712 compute  dask/base.py:603
   │                       [7 frames hidden]  dask, queue, threading, <built-in>
   │                          100.197 lock.acquire  <built-in>
   │        15.078 JointLikelihood.minus_log_like_profile  threeML/classicMLE/joint_likelihood.py:936
   │        └─ 15.078 HAL.inner_fit  hawc_hal/HAL.py:1609
   │           └─ 15.078 HAL.get_log_like  hawc_hal/HAL.py:1019
   │              └─ 15.078 HAL._multithreading_log_like  hawc_hal/HAL.py:974
   │                 └─ 15.078 compute  dask/base.py:603
   │                       [7 frames hidden]  dask, queue, threading, <built-in>
   ├─ 34.655 HAL.__init__  hawc_hal/HAL.py:73
   │  ├─ 20.678 map_tree_factory  hawc_hal/maptree/map_tree.py:21
   │  │  └─ 20.678 MapTree.from_root_file  hawc_hal/maptree/map_tree.py:54
   │  │     └─ 20.678 from_root_file  hawc_hal/maptree/from_root_file.py:165
   │  │        └─ 20.678 Pool.starmap  multiprocessing/pool.py:369
   │  │              [6 frames hidden]  multiprocessing, threading, <built-in>
   │  └─ 13.475 hawc_response_factory  hawc_hal/response/response.py:31
   │     └─ 13.475 HAWCResponse.from_root_file  hawc_hal/response/response.py:381
   │        ├─ 9.500 ResponseBin.from_ttree  hawc_hal/response/response_bin.py:172
   │        │  └─ 7.750 PSFWrapper.psf_eval  hawc_hal/psf_fast/psf_wrapper.py:222
   │        │     ├─ 4.500 PSFWrapper.__init__  hawc_hal/psf_fast/psf_wrapper.py:50
   │        │     │  └─ 3.750 PSFWrapper._prepare_brightness_interpolation_points  hawc_hal/psf_fast/psf_wrapper.py:94
   │        │     │     └─ 3.750 <listcomp>  hawc_hal/psf_fast/psf_wrapper.py:100
   │        │     │        └─ 2.750 PSFWrapper.integral  hawc_hal/psf_fast/psf_wrapper.py:252
   │        │     │           └─ 2.750 InterpolatedUnivariateSpline.integral  scipy/interpolate/_fitpack2.py:423
   │        │     │                 [2 frames hidden]  scipy
   │        │     └─ 2.250 <listcomp>  hawc_hal/psf_fast/psf_wrapper.py:242
   │        │        └─ 2.000 psf_func  hawc_hal/psf_fast/psf_wrapper.py:195
   │        └─ 3.725 Pool.starmap  multiprocessing/pool.py:369
   │              [6 frames hidden]  multiprocessing, threading, <built-in>
   ├─ 13.867 JointLikelihood.compute_TS  threeML/classicMLE/joint_likelihood.py:1289
   │     [2 frames hidden]  threeML
   │        13.616 JointLikelihood._assign_model_to_data  threeML/classicMLE/joint_likelihood.py:118
   │        └─ 13.616 HAL.set_model  hawc_hal/HAL.py:411
   │           └─ 13.616 HAL._setup_psf_convolutors  hawc_hal/HAL.py:257
   │              └─ 13.616 PSFConvolutor.__init__  hawc_hal/psf_fast/psf_convolutor.py:18
   │                 └─ 13.616 PSFInterpolator.point_source_image  hawc_hal/psf_fast/psf_interpolator.py:123
   │                    └─ 13.616 PSFInterpolator._get_point_source_image_aitoff  hawc_hal/psf_fast/psf_interpolator.py:36
   │                       └─ 13.114 reproject_exact  reproject/spherical_intersect/high_level.py:11
   │                             [4 frames hidden]  reproject
   └─ 13.711 JointLikelihood.__init__  threeML/classicMLE/joint_likelihood.py:52
         [2 frames hidden]  threeML
            13.711 JointLikelihood._assign_model_to_data  threeML/classicMLE/joint_likelihood.py:118
            └─ 13.711 HAL.set_model  hawc_hal/HAL.py:411
               └─ 13.711 HAL._setup_psf_convolutors  hawc_hal/HAL.py:257
                  └─ 13.711 PSFConvolutor.__init__  hawc_hal/psf_fast/psf_convolutor.py:18
                     └─ 13.711 PSFInterpolator.point_source_image  hawc_hal/psf_fast/psf_interpolator.py:123
                        └─ 13.711 PSFInterpolator._get_point_source_image_aitoff  hawc_hal/psf_fast/psf_interpolator.py:36
                           └─ 13.209 reproject_exact  reproject/spherical_intersect/high_level.py:11
                                 [4 frames hidden]  reproject