Olga Sergijenko / ctools

Commit abb29aa8bbf3336fcd96b6470f086e97f72095db

Authored by Luigi Tibaldo
1 parent 8ebb4df8

csphasecrv parallelized

Inline Side-by-side
Showing 2 changed files with 177 additions and 76 deletions
cscripts/csphasecrv.par
  View file @abb29aa
... ... @@ -55,12 +55,13 @@ statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
55 55 #
56 56 # Standard parameters
57 57 #====================
58   -publish, b, h, no,,, "Publish phase curve on VO Hub?"
59   -chatter, i, h, 2,0,4, "Chattiness of output"
60   -clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
61   -debug, b, h, no,,, "Activate debugging mode?"
62   -mode, s, h, ql,,, "Mode of automatic parameters"
63   -logfile, f, h, csphasecrv.log,,, "Log filename"
  58 +nthreads, i, h, 0,0,1000, "Number of parallel processes (0=use all available CPUs)"
  59 +publish, b, h, no,,, "Publish phase curve on VO Hub?"
  60 +chatter, i, h, 2,0,4, "Chattiness of output"
  61 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  62 +debug, b, h, no,,, "Activate debugging mode?"
  63 +mode, s, h, ql,,, "Mode of automatic parameters"
  64 +logfile, f, h, csphasecrv.log,,, "Log filename"
64 65  
65 66  
66 67  
... ...
cscripts/csphasecrv.py
  View file @abb29aa
... ... @@ -23,6 +23,33 @@ import gammalib
23 23 import ctools
24 24 from cscripts import obsutils
25 25 from cscripts import ioutils
  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, phbin):
  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._phase_bin(phbin)
  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 # ================ #
... ... @@ -71,6 +98,50 @@ class csphasecrv(ctools.csobservation):
71 98 # Return
72 99 return
73 100  
  101 + # State methods por pickling
  102 + def __getstate__(self):
  103 + """
  104 + Extend ctools.csobservation getstate method to include some members
  105 +
  106 + Returns
  107 + -------
  108 + state : dict
  109 + Pickled instance
  110 + """
  111 + # Set pickled dictionary
  112 + state = {'base' : ctools.csobservation.__getstate__(self),
  113 + 'srcname' : self._srcname,
  114 + 'phbins' : self._phbins,
  115 + 'stacked' : self._stacked,
  116 + 'onoff' : self._onoff,
  117 + 'fits' : self._fits,
  118 + 'fitmodels': self._fitmodels,
  119 + 'nthreads' : self._nthreads}
  120 +
  121 + # Return pickled dictionary
  122 + return state
  123 +
  124 + def __setstate__(self, state):
  125 + """
  126 + Extend ctools.csobservation setstate method to include some members
  127 +
  128 + Parameters
  129 + ----------
  130 + state : dict
  131 + Pickled instance
  132 + """
  133 + ctools.csobservation.__setstate__(self, state['base'])
  134 + self._srcname = state['srcname']
  135 + self._phbins = state['phbins']
  136 + self._onoff = state['onoff']
  137 + self._stacked = state['stacked']
  138 + self._fits = state['fits']
  139 + self._fitmodels= state['fitmodels']
  140 + self._nthreads = state['nthreads']
  141 +
  142 + # Return
  143 + return
  144 +
74 145 # Private methods
75 146 def _get_parameters(self):
76 147 """
... ... @@ -111,6 +182,9 @@ class csphasecrv(ctools.csobservation):
111 182 # Write input parameters into logger
112 183 self._log_parameters(gammalib.TERSE)
113 184  
  185 + # Set number of processes for multiprocessing
  186 + self._nthreads = mputils.nthreads(self)
  187 +
114 188 # Return
115 189 return
116 190  
... ... @@ -320,6 +394,79 @@ class csphasecrv(ctools.csobservation):
320 394 # Return
321 395 return
322 396  
  397 + def _phase_bin(self,phbin):
  398 +
  399 + # Write time bin into header
  400 + self._log_header2(gammalib.TERSE, 'PHASE %f - %f' %
  401 + (phbin[0], phbin[1]))
  402 +
  403 + # Select events
  404 + select = ctools.ctselect(self.obs().copy())
  405 + select['emin'] = self['emin'].real()
  406 + select['emax'] = self['emax'].real()
  407 + select['tmin'] = 'UNDEFINED'
  408 + select['tmax'] = 'UNDEFINED'
  409 + select['rad'] = 'UNDEFINED'
  410 + select['ra'] = 'UNDEFINED'
  411 + select['dec'] = 'UNDEFINED'
  412 + select['expr'] = 'PHASE>' + str(phbin[0]) + ' && PHASE<' + str(phbin[1])
  413 + select.run()
  414 +
  415 + # Set phase string
  416 + phstr = str(phbin[0]) + '-' + str(phbin[1])
  417 +
  418 + # Add phase to observation id
  419 + for i in range(0, select.obs().size()):
  420 + oldid = select.obs()[i].id()
  421 + select.obs()[i].id(oldid + '_' + phstr)
  422 + obs = select.obs()
  423 +
  424 + # If an On/Off analysis is requested generate the On/Off observations
  425 + if self._onoff:
  426 + obs = obsutils.get_onoff_obs(self, select.obs())
  427 +
  428 + # ... otherwise, if stacked analysis is requested then bin the
  429 + # events and compute the stacked response functions and setup
  430 + # an observation container with a single stacked observation.
  431 + elif self._stacked:
  432 + obs = obsutils.get_stacked_obs(self, select.obs())
  433 +
  434 + # Header
  435 + self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
  436 +
  437 + # The On/Off analysis can produce empty observation containers,
  438 + # e.g., when on-axis observations are used. To avoid ctlike asking
  439 + # for a new observation container (or hang, if in interactive mode)
  440 + # we'll run ctlike only if the size is >0
  441 + if obs.size() > 0:
  442 +
  443 + # Do maximum likelihood model fitting
  444 + like = ctools.ctlike(obs)
  445 + like['edisp'] = self['edisp'].boolean()
  446 + like.run()
  447 +
  448 + # Renormalize models to phase selection
  449 + # TODO move the scaling from the temporal to the spectral component
  450 + for model in like.obs().models():
  451 + scaled_norm = model['Normalization'].value() / (phbin[1] - phbin[0])
  452 + model['Normalization'].value(scaled_norm)
  453 +
  454 + # Store fit model
  455 + fitmodels = like.obs().models().copy()
  456 +
  457 + # ... otherwise we set an empty model container
  458 + else:
  459 + self._log_string(gammalib.TERSE,
  460 + 'PHASE %f - %f: no observations available'
  461 + ' for fitting' % (phbin[0], phbin[1]))
  462 +
  463 + # Set empty models container
  464 + fitmodels = gammalib.GModels()
  465 +
  466 + result = {'phstr' :phstr,
  467 + 'fitmodels' :fitmodels}
  468 +
  469 + return result
323 470  
324 471 # Public methods
325 472 def run(self):
... ... @@ -336,84 +483,37 @@ class csphasecrv(ctools.csobservation):
336 483 # Write observation into logger
337 484 self._log_observations(gammalib.NORMAL, self.obs(), 'Observation')
338 485  
339   - # Store and clear observations to split them into phase bins
340   - orig_obs = self.obs().copy()
341   -
342 486 # Dictionary to save phase fitted models
343 487 self._fitmodels = {}
344 488  
345 489 # Write header
346 490 self._log_header1(gammalib.TERSE, 'Generate phase curve')
347 491  
348   - # Loop over all phases
349   - for phbin in self._phbins:
  492 + # If using multiprocessing
  493 + if self._nthreads > 1:
350 494  
351   - # Write time bin into header
352   - self._log_header2(gammalib.TERSE, 'PHASE %f - %f' %
353   - (phbin[0], phbin[1]))
354   -
355   - # Select events
356   - select = ctools.ctselect(orig_obs)
357   - select['emin'] = self['emin'].real()
358   - select['emax'] = self['emax'].real()
359   - select['tmin'] = 'UNDEFINED'
360   - select['tmax'] = 'UNDEFINED'
361   - select['rad'] = 'UNDEFINED'
362   - select['ra'] = 'UNDEFINED'
363   - select['dec'] = 'UNDEFINED'
364   - select['expr'] = 'PHASE>'+str(phbin[0])+' && PHASE<'+str(phbin[1])
365   - select.run()
366   -
367   - # Set phase string
368   - phstr = str(phbin[0]) + '-' + str(phbin[1])
369   -
370   - # Add phase to observation id
371   - for i in range(0,select.obs().size()):
372   - oldid = select.obs()[i].id()
373   - select.obs()[i].id(oldid+'_'+phstr)
374   - obs = select.obs()
375   -
376   - # If an On/Off analysis is requested generate the On/Off observations
377   - if self._onoff:
378   - obs = obsutils.get_onoff_obs(self, select.obs())
379   -
380   - # ... otherwise, if stacked analysis is requested then bin the
381   - # events and compute the stacked response functions and setup
382   - # an observation container with a single stacked observation.
383   - elif self._stacked:
384   - obs = obsutils.get_stacked_obs(self, select.obs())
385   -
386   - # Header
387   - self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
388   -
389   - # The On/Off analysis can produce empty observation containers,
390   - # e.g., when on-axis observations are used. To avoid ctlike asking
391   - # for a new observation container (or hang, if in interactive mode)
392   - # we'll run ctlike only if the size is >0
393   - if obs.size() > 0:
394   -
395   - # Do maximum likelihood model fitting
396   - like = ctools.ctlike(obs)
397   - like['edisp'] = self['edisp'].boolean()
398   - like.run()
399   -
400   - # Renormalize models to phase selection
401   - # TODO move the scaling from the temporal to the spectral component
402   - for model in like.obs().models():
403   - scaled_norm = model['Normalization'].value()/(phbin[1]-phbin[0])
404   - model['Normalization'].value(scaled_norm)
405   -
406   - # Store fit model
407   - self._fitmodels[phstr] = like.obs().models().copy()
408   -
409   - # ... otherwise we set an empty model container
410   - else:
411   - self._log_string(gammalib.TERSE,
412   - 'PHASE %f - %f: no observations available'
413   - ' for fitting' %(phbin[0], phbin[1]))
  495 + # Create pool of workers
  496 + from multiprocessing import Pool
  497 + pool = Pool(processes = self._nthreads)
  498 +
  499 + # Run time bin analysis in parallel with map
  500 + args = [(self, phbin) for phbin in self._phbins]
  501 + poolresults = pool.map(_multiprocessing_func_wrapper, args)
  502 +
  503 + # Close pool and join
  504 + pool.close()
  505 + pool.join()
  506 +
  507 + # Construct results
  508 + for i in range(len(self._phbins)):
  509 + self._fitmodels[poolresults[i][0]['phstr']] = poolresults[i][0]['fitmodels']
  510 + self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False)
414 511  
415   - # Set empty models container
416   - self._fitmodels[phstr] = gammalib.GModels()
  512 + # Otherwise, loop over all phases
  513 + else:
  514 + for phbin in self._phbins:
  515 + result = self._phase_bin(phbin)
  516 + self._fitmodels[result['phstr']] = result['fitmodels']
417 517  
418 518 # Create FITS file
419 519 self._create_fits()
... ...