Iftach Sadeh / ctools

Commit 2a1dd3adc80092132c440d05010157bbfdcd03bd

Authored by Luigi Tibaldo
1 parent 56569b2e

parallelisation

Inline Side-by-side
Showing 2 changed files with 186 additions and 33 deletions
cscripts/csphagen.par
  View file @2a1dd3a
... ... @@ -53,8 +53,9 @@ etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
53 53 #
54 54 # Standard parameters
55 55 #====================
56   -chatter, i, h, 2,0,4, "Chattiness of output"
57   -clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
58   -debug, b, h, no,,, "Activate debugging mode?"
59   -mode, s, h, ql,,, "Mode of automatic parameters"
60   -logfile, f, h, csphagen.log,,, "Log filename"
  56 +nthreads, i, h, 0,0,1000, "Number of parallel processes (0=use all available CPUs)"
  57 +chatter, i, h, 2,0,4, "Chattiness of output"
  58 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  59 +debug, b, h, no,,, "Activate debugging mode?"
  60 +mode, s, h, ql,,, "Mode of automatic parameters"
  61 +logfile, f, h, csphagen.log,,, "Log filename"
... ...
cscripts/csphagen.py
  View file @2a1dd3a
... ... @@ -23,6 +23,33 @@ import gammalib
23 23 import ctools
24 24 import math
25 25 import sys
  26 +from cscripts import mputils
  27 +
  28 +# ============================================ #
  29 +# Global functions for multiprocessing support #
  30 +# ============================================ #
  31 +def _multiprocessing_func_wrapper(args):
  32 + return _multiprocessing_func(*args)
  33 +def _multiprocessing_func(cls, i):
  34 +
  35 + # Initialise thread logger
  36 + cls._log.clear()
  37 + cls._log.buffer_size(100000)
  38 +
  39 + # Compute light curve bin
  40 + cstart = cls.celapse()
  41 + result = cls._process_observation(i)
  42 + celapse = cls.celapse() - cstart
  43 + buffer = cls._log.buffer()
  44 +
  45 + # Close logger
  46 + cls._log.close()
  47 +
  48 + # Collect thread information
  49 + info = {'celapse': celapse, 'log': buffer}
  50 +
  51 + # Return light curve bin result and thread information
  52 + return result, info
26 53  
27 54  
28 55 # =============== #
... ... @@ -43,6 +70,7 @@ class csphagen(ctools.csobservation):
43 70  
44 71 # Initialise other variables
45 72 self._ebounds = gammalib.GEbounds()
  73 + self._etruebounds = gammalib.GEbounds()
46 74 self._src_dir = gammalib.GSkyDir()
47 75 self._src_reg = gammalib.GSkyRegions()
48 76 self._src_model = None
... ... @@ -51,10 +79,62 @@ class csphagen(ctools.csobservation):
51 79 self._has_exclusion = False
52 80 self._srcshape = ''
53 81 self._rad = 0.0
  82 + self._nthreads = 0
54 83  
55 84 # Return
56 85 return
57 86  
  87 + # State methods por pickling
  88 + def __getstate__(self):
  89 + """
  90 + Extend ctools.csobservation getstate method to include some members
  91 +
  92 + Returns
  93 + -------
  94 + state : dict
  95 + Pickled instance
  96 + """
  97 + # Set pickled dictionary
  98 + state = {'base' : ctools.csobservation.__getstate__(self),
  99 + 'ebounds' : self._ebounds,
  100 + 'etruebounds' : self._etruebounds,
  101 + 'src_dir' : self._src_dir,
  102 + 'src_reg' : self._src_reg,
  103 + 'src_model' : self._src_model,
  104 + 'bkg_regs' : self._bkg_regs,
  105 + 'excl_reg' : self._excl_reg,
  106 + 'has_exclusion': self._has_exclusion,
  107 + 'srcshape' : self._srcshape,
  108 + 'rad' : self._rad,
  109 + 'nthreads' : self._nthreads}
  110 +
  111 + # Return pickled dictionary
  112 + return state
  113 +
  114 + def __setstate__(self, state):
  115 + """
  116 + Extend ctools.csobservation setstate method to include some members
  117 +
  118 + Parameters
  119 + ----------
  120 + state : dict
  121 + Pickled instance
  122 + """
  123 + ctools.csobservation.__setstate__(self, state['base'])
  124 + self._ebounds = state['ebounds']
  125 + self._etruebounds = state['etruebounds']
  126 + self._src_dir = state['src_dir']
  127 + self._src_reg = state['src_reg']
  128 + self._src_model = state['src_model']
  129 + self._bkg_regs = state['bkg_regs']
  130 + self._excl_reg = state['excl_reg']
  131 + self._has_exclusion = state['has_exclusion']
  132 + self._srcshape = state['srcshape']
  133 + self._rad = state['rad']
  134 + self._nthreads = state['nthreads']
  135 +
  136 + # Return
  137 + return
58 138  
59 139 # Private methods
60 140 def _query_src_direction(self):
... ... @@ -205,6 +285,9 @@ class csphagen(ctools.csobservation):
205 285 # Write input parameters into logger
206 286 self._log_parameters(gammalib.TERSE)
207 287  
  288 + # Set number of processes for multiprocessing
  289 + self._nthreads = mputils.nthreads(self)
  290 +
208 291 # Return
209 292 return
210 293  
... ... @@ -378,6 +461,80 @@ class csphagen(ctools.csobservation):
378 461 # Return background regions
379 462 return bkg_reg
380 463  
  464 + def _process_observation(self,i):
  465 + """
  466 + Generate On/Off spectra for individual observation
  467 +
  468 + :param i: int, observation number
  469 + :return: dictionary
  470 + On/Off spectra, background regions, observation id
  471 + """
  472 +
  473 + # Retrieve observation from container
  474 + obs = self.obs()[i]
  475 +
  476 + # Skip non CTA observations
  477 + if obs.classname() != 'GCTAObservation':
  478 + self._log_string(gammalib.NORMAL, 'Skip %s observation "%s"' % \
  479 + (obs.instrument(), obs.id()))
  480 + onoff = 0
  481 + bkg_reg = 0
  482 +
  483 + # Otherwise calculate On/Off spectra
  484 + else:
  485 + # Log current observation
  486 + self._log_string(gammalib.NORMAL, 'Process observation "%s"' % \
  487 + (obs.id()))
  488 +
  489 + # Set background regions for this observation
  490 + bkg_reg = self._set_background_regions(obs)
  491 +
  492 + # If there are background regions then create On/Off observation
  493 + # and append it to the output container
  494 + if bkg_reg.size() >= self['bkgregmin'].integer():
  495 + onoff = gammalib.GCTAOnOffObservation(obs,
  496 + self._src_model,
  497 + self._etruebounds,
  498 + self._ebounds,
  499 + self._src_reg,
  500 + bkg_reg)
  501 + onoff.id(obs.id())
  502 +
  503 + else:
  504 + self._log_string(gammalib.NORMAL, 'Observation %s not included '
  505 + 'in spectra generation' % (obs.id()))
  506 +
  507 + # Construct dictionary with results
  508 + result = {'onoff' : onoff,
  509 + 'bkg_reg' : bkg_reg,
  510 + 'id' : obs.id()}
  511 +
  512 + # Return results
  513 + return result
  514 +
  515 + def _unpack_result(self,outobs,result):
  516 + """
  517 + Unpack result from calculation of On/Off regions and appends it to observation
  518 + container.
  519 +
  520 + :param outobs: `~gammalib.GObservations`, observation container
  521 + :param result: dict, result to unpack
  522 + :return: `~gammalib.GObservations`
  523 + observation container with result appended
  524 + """
  525 +
  526 + # If the results contain an On/Off observation
  527 + if result['onoff'].classname() == 'GCTAOnOffObservation':
  528 + # Append observation to observation container
  529 + outobs.append(result['onoff'])
  530 +
  531 + # Append background regions
  532 + self._bkg_regs.append({'regions': result['bkg_reg'],
  533 + 'id': result['id']})
  534 +
  535 + # Return observation container
  536 + return outobs
  537 +
381 538  
382 539 # Public methods
383 540 def run(self):
... ... @@ -395,7 +552,7 @@ class csphagen(ctools.csobservation):
395 552 self._log_observations(gammalib.NORMAL, self.obs(), 'Observation')
396 553  
397 554 # Set true energy bins
398   - etrue_ebounds = self._etrue_ebounds()
  555 + self._etruebounds = self._etrue_ebounds()
399 556  
400 557 # Write header
401 558 self._log_header1(gammalib.TERSE, 'Spectral binning')
... ... @@ -418,37 +575,32 @@ class csphagen(ctools.csobservation):
418 575 if self._src_model == None:
419 576 self._src_model = gammalib.GModelSpatialPointSource(self._src_dir)
420 577  
421   - # Loop through observations and generate pha, arf, rmf files
422   - for obs in self.obs():
  578 + # If there is more than one observation and we use multiprocessing
  579 + if self._nthreads > 1 and self.obs().size() > 1:
423 580  
424   - # Skip non CTA observations
425   - if obs.classname() != 'GCTAObservation':
426   - self._log_string(gammalib.NORMAL, 'Skip %s observation "%s"' % \
427   - (obs.instrument(), obs.id()))
428   - continue
  581 + # Create pool of workers
  582 + from multiprocessing import Pool
  583 + pool = Pool(processes = self._nthreads)
429 584  
430   - # Log current observation
431   - self._log_string(gammalib.NORMAL, 'Process observation "%s"' % \
432   - (obs.id()))
  585 + # Run time bin analysis in parallel with map
  586 + args = [(self, i) for i in range(self.obs().size())]
  587 + poolresults = pool.map(_multiprocessing_func_wrapper, args)
433 588  
434   - # Set background regions for this observation
435   - bkg_reg = self._set_background_regions(obs)
  589 + # Close pool and join
  590 + pool.close()
  591 + pool.join()
436 592  
437   - # If there are background regions then create On/Off observation
438   - # and append it to the output container
439   - if bkg_reg.size() >= self['bkgregmin'].integer():
440   - onoff = gammalib.GCTAOnOffObservation(obs,
441   - self._src_model,
442   - etrue_ebounds,
443   - self._ebounds,
444   - self._src_reg,
445   - bkg_reg)
446   - onoff.id(obs.id())
447   - outobs.append(onoff)
448   - self._bkg_regs.append({'regions': bkg_reg, 'id': obs.id()})
449   - else:
450   - self._log_string(gammalib.NORMAL, 'Observation %s not included '
451   - 'in spectra generation' % (obs.id()))
  593 + # Construct results
  594 + for i in range(self.obs().size()):
  595 + outobs = self._unpack_result(outobs,poolresults[i][0])
  596 + self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False)
  597 +
  598 + # Otherwise, loop through observations and generate pha, arf, rmf files
  599 + else:
  600 + for i in range(self.obs().size()):
  601 + # Process individual observation
  602 + result = self._process_observation(i)
  603 + outobs = self._unpack_result(outobs, result)
452 604  
453 605 # Stack observations
454 606 if outobs.size() > 1 and self['stack'].boolean() == True:
... ...