Source code for tesfdmtools.methods.aldtimes

#!/usr/bin/env python3
'''
.. module:: alltimes
   :synopsis: Retrieve delta time to previous event considering events of all pixels  
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>

'''
import numpy

from mpi4py import MPI
from tesfdmtools.utils.cu import cu

def aldtimes(comm,rank,hdf,channel=None,frequencies=None,nrecs="all",\
                threshold=2000,debug=False):
   '''
      Go through all frequencies and store eventtimes. Combine and sort these times to
      obtain delta times for all events to the previous event, and separate for
      each frequency (pixel).
      
      Args:
         * `comm` = the MPI identifier for the common process
         * `rank` = the rank of this process
         *  `hdf` = the hdf data file
        
      Kwargs:
         *     `channel` = the hdf channel to use
         * `frequencies` = list of frequencies to process for this pid
         *       `nrecs` = number of records to process ( **all** is all records )
         *   `threshold` = threshold for events to consider
         *       `debug` = debug flag; when **True** show debug information 
   '''
   
   if channel is None:
      channel=hdf.channels[0]
            
   if frequencies is None:
      frequencies=[hdf.channel[channel].freq_conf['pixel_index'][0]]
      
   evtimes={}
   trecs=0
         
   for j,freq in frequencies:
         
      print("aldtimes: read records from channel %d, frequency number %d" % (channel,freq))

      nrecords=cu(hdf.channel[channel].attrs['Nevents'])

      if ( type(nrecs) is str ) and ( nrecs == 'all' ):
         nrecs=nrecords
      else:
         nrecs=min(nrecords,int(nrecs))
      trecs=trecs+nrecs
      mrecs=nrecs-1
      
      evtimes[freq]=numpy.zeros(nrecs,dtype=float)
      sirecord=hdf.channel[channel].freq[freq][0][:,0]
      samplefreq=float(hdf.channel[channel].freq[freq].attrs['sample_rate'])  # sample frequency
 
      for i,record in enumerate(hdf.channel[channel].freq[freq]):

         irecord=record[:,0]                                # use I-signal only
         
         ptp=numpy.ptp(irecord)
         if ptp > threshold:
            rectim=float(record.attrs['event_time'])        # time of record
            position=numpy.argmin(irecord)
            evtimes[freq][i]=rectim+position/samplefreq # time of the event position in seconds
            
         if i >= mrecs:
            break
            
      if debug:
         np=20
         ff=("freq: %d times:" % freq)+np*" %6.2f"+"\n"
         print (ff % tuple(evtimes[freq][0:np]) )
         
            
   aevtimes=comm.gather([trecs,evtimes],root=0) # get all eventtimes from all processes together
   
   if rank == 0:
   
      findxs={}					# keep indices of frequencies for grand time table
      tabindx=0
      
      js=[]					# process numbers for all frequencies
      ifreq=[]					# all frequencies
      
      trec=0
      for j,drec in enumerate(aevtimes):
         trec=trec+drec[0]      		# total number of events
         if debug:
            print("Number of events for this process: ",drec[0])
         for freq in drec[1]:
            js.append(j)
            ifreq.append(freq)
            
      sjindx=numpy.argsort(ifreq)   # sort on frequency, to obtain identical
                                    # frequency order, independant of process order
                                    # otherwise, sorting of the times may yield slightly different results
         
      if debug:
         print("Total number of events: ",trec)
            
      etimesl=numpy.zeros(trec,dtype=float)	# grant total of all event times
      for sj in sjindx:
         drec=aevtimes[js[sj]]
         freq=ifreq[sj]
         nrecs=drec[1][freq].size		# number of events for this frequency
         if debug:
            print("number of events for frequency: ",freq,"  n: ",nrecs)
         findxs[freq]=[tabindx,tabindx+nrecs] # hold indices in list for this frequency
         etimesl[tabindx:(tabindx+nrecs)]=drec[1][freq] # add event times to list
         tabindx=tabindx+nrecs
            
      tindx=numpy.argsort(etimesl)		# sort on time order
      pindx=numpy.concatenate(([tindx[0]],tindx[:-1]))
      deltap=numpy.zeros(tindx.size)
      deltap[tindx]=etimesl[tindx]-etimesl[pindx] # delta times to previous event
      
      if debug:
         np=20
         ff=np*" %6.2f"
         print (("alltimes:"+ff+"\n") % tuple(etimesl[tindx[0:np]]) )
         print (("alldelta:"+ff+"\n") % tuple(deltap[tindx][0:np]) )
      
      paltimes=[]				# assemble data for scatter to the processes
      for drec in aevtimes:
         pevtimes={}
         for freq in drec[1]:
            pevtimes[freq]=deltap[findxs[freq][0]:findxs[freq][1]]
            
            if debug:
               np=20
               ff=("freq: %d deltas" % freq)+np*" %6.2f"
               print (ff % tuple(pevtimes[freq][0:np]) ) 
               print ("freq: %d ndeltas: %d" % (freq,pevtimes[freq].size) )
            
         paltimes.append(pevtimes)
         
   else:
   
      paltimes=None
      
   pevtimes=comm.scatter(paltimes,root=0)
   
   return pevtimes
   
#===================================================================

def selectald(indx,deltatims,etprev=[None,None],debug=False):
   '''
      Make new index for events which have delta-times within the given window
      
      Args:
         *      `indx` = index of events to consider
         * `deltatims` = delta times to previous event for all events
         
      Kwargs:
         * `etprev` = time window [min,max] for valid event
         *  `debug` = print debug information when set **True**
   '''

   if etprev[0] is None:
      return indx
   
   if debug:
      np=20
      pfmt="selectald: deltas: "+np*" %6.2f"+"\n"
      print(pfmt % tuple(deltatims[0:20]) )
      
   sindx,=numpy.where( ( deltatims[indx] > etprev[0] ) & \
                       ( deltatims[indx] < etprev[1] ) )
                       
#   print("selectald: n(indx) n(deltatims) n(sindx): ",indx.size,deltatims.size,sindx.size)
                       
   nindx=indx[sindx]
   
   return nindx