;+
; NAME:
;   GTIMERGE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Merge two Good Time Interval (GTIs) arrays into a single array
;
; CALLING SEQUENCE:
;   NEWGTI = GTIMERGE(GTI1, GTI2, COUNT=COUNT, [/INTERSECT, /UNION,
;                     /INVERT1, /INVERT2, TTOLERANCE=])
;
; DESCRIPTION: 
;
;   The function GTIMERGE accepts two existing valid Good Time
;   Interval (GTI) arrays and merges them into a single array.  Either
;   the intersection or the union of the two GTIs are returned.
;
;   The intersection refers to the set of intervals which lie in both
;   intervals.  The union refers to the set of intervals which lie in
;   at least one but not necessarily both.  Here is an example of both
;   kinds of operations.  Let us start with two intervals here:
;
;        0              50              100                  170   GTI1
;   <----|==============|----------------|====================|------>
;
;                 30                         120                   GTI2
;   <--------------|==========================|---------------------->
;
;   These intervals would be represented by GTI1=[[0,50],[100,170]]
;   and GTI2=[[30,120]]. The intersection of the two sets of intervals
;   are the points which lie in both, ie [[30,50],[100,120]]:
;
;                 30   50               100  120              INTERSECT
;   <--------------|====|----------------|====|---------------------->
;
;   The union is the combination of both intervals, ie [[0,170]]:
;
;        0                                                   170  UNION
;   <----|====================================================|------>
;
;   It is also possible to treat either one of the input arrays as
;   "bad" intervals using the INVERT1 and/or INVERT2 keywords.  When
;   an interval is inverted, then the output is composed only of areas
;   *outside* the specified intervals.
;
;   It should be noted that this function is not constrained to
;   operation only on time arrays.  It should work on any
;   one-dimensional quantity with intervals.
;
;
; PERFORMANCE:  Combining many intervals
;
;   Users who wish to combine many intervals in sequence will find a
;   performance degradation.  The problem is that each GTIMERGE
;   operation is order N^2 execution time where N is the number of
;   intervals.  Thus, if N mostly distinct GTIs are merged, then the
;   running time will be order N^3.  This is unacceptable, but there
;   is a workaround.
;
;   Users can accumulate "sub" GTIs by merging subsets of the full
;   number of intervals to be merged, and then occasionally merging
;   into the final output GTI.  As an example, here first is a simple
;   merging of 1000 different GTIs:
;
;     totgti = 0L                   ;; Empty GTI
;     FOR i = 0, 999 DO BEGIN
;         gti = ...                
;         totgti = gtimerge(totgti, gti, /union)
;     ENDFOR
;
;   This computation may take a long time.  Instead the merging can be
;   broken into chunks.
;
;     totgti = 0L
;     chgti  = 0L      ;; "Chunk" GTI
;     FOR i = 0, 999 DO BEGIN
;         gti = ...
;         chgti = gtimerge(chgti, gti, /union)
;         if (n_elements(chgti) GT 100) OR (i EQ 999) then begin
;             ;; Merge "chunk" gti into final one, and reset
;             totgti = gtimerge(totgti, chgti, /union)
;             chgti = 0L
;         endif
;     ENDFOR
;
;   Note that the final merge is guaranteed because of the (i EQ 999)
;   comparison.
;
; INPUTS:
;
;   GTI1, GTI2 - the two input GTI arrays.
;
;         Each array is a 2xNINTERVAL array where NINTERVAL is the
;         number of intervals, which can be different for each array.
;         GTI(*,i) represents the start and stop times of interval
;         number i.  The intervals must be non-overlapping and
;         time-ordered (use GTITRIM to achieve this).
;
;         A scalar value of zero indicates that the GTI is empty, ie,
;         there are no good intervals.
;
; KEYWORDS:
;
;   INTERSECT - if set, then the resulting GTI contains only those
;               intervals that are in both input sets.
;
;   UNION - if set, then the resulting GTI contains those intervals
;           that are in either input set.
;
;   COUNT - upon return, the number of resulting intervals.  A value
;           of zero indicates no good time intervals.
;
;   INVERT1 - if set, then GTI1 is considered to be inverted, ie, a
;             set of "bad" intervals rather than good.
;
;   INVERT2 - if set, then GTI2 is considered to be inverted, ie, a
;             set of "bad" intervals rather than good.
;
;   TTOLERANCE - a scalar value indicating the tolerance for
;                determining whether values are equal.  This number
;                can be important for intervals that do not match
;                precisely.
;                Default: Machine precision
;
; RETURNS:
;
;   A new GTI array containing the merged intervals.  The array is
;   2xCOUNT where COUNT is the number of resulting intervals.
;   GTI(*,i) represents the start and stop times of interval number i.
;   The intervals are non-overlapping and time-ordered.
;
;   If COUNT is zero then the returned array is a scalar value of
;   zero, indicating no good intervals were found.
;
; SEE ALSO:
;
;   GTITRIM, GTIENLARGE
;
; MODIFICATION HISTORY:
;   Written, CM, 1997-2001
;   Documented, CM, Apr 2001
;   Handle case of zero-time GTIs, CM, 02 Aug 2001
;   Handle "fractured" GTIs correctly, though worriedly, CM, 15 Oct
;     2001
;   Handle case where both inputs are empty, but /INVERT1 and/or
;     /INVERT2 are set, CM, 08 Aug 2006
;
;  $Id: gtimerge.pro,v 1.7 2006/10/22 09:49:53 craigm Exp $
;
;-
; Copyright (C) 1997-2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function gtimerge, gti1, gti2, count=count, intersect=intersect, union=union, $
                   invert1=invert1, invert2=invert2, quiet=quiet, $
                   ttolerance=ttol, query=query

  if keyword_set(query) then return, 1

  ;; Intialize count variable.  None found yet.
  count = 0L

  forward_function gtitrim
  catch, catcherr & qgti = 0
  if catcherr EQ 0 then qgti = gtitrim(/query)
  catch, /cancel
  if catcherr NE 0 OR qgti EQ 0 then $
    message, 'ERROR: The function GTITRIM must be in your IDL path'

  if n_elements(ttol) EQ 0 then begin
      sz1 = size(gti1) & sz2 = size(gti2)
      if sz1(sz1(0)+1) EQ 5 OR sz2(sz2(0)+1) EQ 5 then $
        double = 1 $
      else $
        double = 0
      mach = machar(double=double)
      ttol = mach.eps
  endif
  if keyword_set(intersect) and keyword_set(union) then $
    message, 'ERROR: cannot perform both INTERSECT and UNION.'

  if NOT keyword_set(intersect) and NOT keyword_set(union) then $
    message, 'ERROR: must perform either INTERSECT or UNION.'

  ;; The strategy is as follows: break the two GTIs down into a list
  ;; of times.  Associated with each time is the ON/OFF status of the
  ;; GTI.  Actually, it's ON, OFF, or NO CHANGE, since one GTI may
  ;; change when the other does not.  The list of times is sorted, and
  ;; then walked through one at a time, in order of time.  When the
  ;; merging conditions are satisfied (either intersect or union),
  ;; then a new gti entry is created.

  count = 0L
  ngti1 = 0L 
  ngti2 = 0L
  if n_elements(gti1) GT 1 then begin
      newgti1 = reform(gti1, 2, n_elements(gti1)/2)
      ;; Normalize and repack the GTI.  First order of business is to
      ;; remove any zero-time intervals.
      newgti1 = gtitrim(newgti1, mingti=ttol, maxgap=ttol, count=ngti1)

      ;; This is the ON/OFF array for gti1
      if ngti1 GT 0 then begin
          on1 = newgti1
          on1(0,*) = 1 - keyword_set(invert1)
          on1(1,*) = 0 + keyword_set(invert1)
      endif
  endif  
  if n_elements(gti2) GT 1 then begin
      newgti2 = reform(gti2, 2, n_elements(gti2)/2)
      ;; Normalize and repack the GTI.  First order of business is to
      ;; remove any zero-time intervals.
      newgti2 = gtitrim(newgti2, mingti=ttol, maxgap=ttol, count=ngti2)

      ;; This is the ON/OFF array for gti1
      if ngti2 GT 0 then begin
          on2 = newgti2
          on2(0,*) = 1 - keyword_set(invert2)
          on2(1,*) = 0 + keyword_set(invert2)
      endif
  endif  

  ;; Merge the one or two ON/OFF arrays together
  if ngti1 GT 0 AND ngti2 GT 0 then begin
      ;; Here is the merged list of times, eventually called tt
      tt = [newgti1(*), newgti2(*)]
      nt = n_elements(tt)
      ;; Now the ON/OFF states are merged as well.  Here, 0 and 1
      ;; represent off and on, and -1 represents NO CHANGE. 
      o1 = reform([on1(*), on2(*)*0-1], nt)
      o2 = reform([on1(*)*0-1, on2(*)], nt)

  endif else if ngti1 GT 0 then begin
      tt = [newgti1(*), min(newgti1)]
      nt = n_elements(tt)
      o1 = reform([on1(*),    -1                   ], nt)
      o2 = reform([on1(*)*0-1, keyword_set(invert2)], nt)
      
  endif else if ngti2 GT 0 then begin
      tt = [min(newgti2), newgti2(*)]
      nt = n_elements(tt)
      o1 = reform([keyword_set(invert1), on2(*)*0-1], nt)
      o2 = reform([-1,                   on2(*)    ], nt)

  endif else begin
      if ((keyword_set(intersect) AND $
           (keyword_set(invert1) AND keyword_set(invert2)))) OR $
        ((keyword_set(union) AND $
          (keyword_set(invert1) OR keyword_set(invert2)))) then begin
          count = 1L
          return, [-!values.f_infinity, +!values.f_infinity]
      endif
      goto, EMPTY_GTI
  endelse
  tt = reform(tt, nt, /overwrite)

  ;; Sort by time
  ts = sort(tt)
  tt = tt(ts) & o1 = o1(ts) & o2 = o2(ts)
  ;; Simulated UNIQ, with well defined properties
  un = where(tt NE shift(tt, +1), ct)
  if ct EQ 0 then wh = 0L
  uu = [un, nt]
  nu = n_elements(un)

  ;; Now we are ready to step through the list of times.  Yes1 and
  ;; yes2 are 1 when the corresponding GTI is on.
  yes1 = keyword_set(invert1) & k1 = NOT keyword_set(invert1)
  yes2 = keyword_set(invert2) & k2 = NOT keyword_set(invert2)
  ;; Within is set to 1 when we are in the middle of a new GTI entry
  within = 0
  numgti = 0L
  dbl1 =  0 & dbl2 =  0
  for j = 0L, nu-1 do begin

      oo1  = -1 &  oo2 = -1
      ;; Collect all times that are the same into one batch.  This is
      ;; important, since a single time point should be considered
      ;; atomic.  If multiple transitions occur in one channel at the
      ;; same time, which shouldn't happen because of the GTI trimming
      ;; above, then we take the first one.
      i = uu(j)
      wh = where(o1(uu(j):uu(j+1)-1) GE 0, ct)
      if ct GT 0 then begin
          o1u = o1(uu(j)+wh)
          oo1 = max(o1u)
          if ct GT 1 then if min(o1u) NE oo1 then dbl1 = 1
      endif
      wh = where(o2(uu(j):uu(j+1)-1) GE 0, ct)
      if ct GT 0 then begin
          o2u = o2(uu(j)+wh)
          oo2 = max(o2u)
          if ct GT 1 then if min(o2u) NE oo2 then dbl2 = 1
      endif
;      print, '  ', o1(uu(j):uu(j+1)-1), format='(A0,100(D4.0," ",:),$)'
;      print, ' --> ', oo1
;      print, '  ', o2(uu(j):uu(j+1)-1), format='(A0,100(D4.0," ",:),$)'
;      print, ' --> ', oo2

      ;; yes1 and yes2 are updated here, if there is a change in gti state
      if oo1 GE 0 then yes1 = (oo1 EQ 1)
      if oo2 GE 0 then yes2 = (oo2 EQ 1)

      ;; Here are the conditions for starting a new GTI entry
      if (NOT within AND (yes1 OR  yes2) AND keyword_set(union)) OR $
        ( NOT within AND (yes1 AND yes2) AND keyword_set(intersect)) then $
        begin
          tstart = tt(i)
          within = 1
      endif
      if dbl1 AND yes1 then yes1 = 0
      if dbl2 AND yes2 then yes1 = 0

      ;; And the conditions for ending it.
      if (within AND (NOT yes1 AND NOT yes2) AND keyword_set(union)) OR $
        ( within AND (NOT yes1 OR  NOT yes2) AND keyword_set(intersect)) then $
        begin
          if numgti EQ 0 then $
            newgti = [ tstart, tt(i) ] $
          else $
            newgti = [ newgti, tstart, tt(i) ]
          numgti   = numgti + 1
          within   = 0

          if dbl1 OR dbl2 then begin
              ;; Now restart the GTI if desired
              tstart = tt(i)
              within = 1
              dbl1 = 0 & dbl2 = 0
          endif
      endif

      LOOPCONT:
  endfor

  ;; I think this clause is activated only if the INVERT1 or INVERT2
  ;; keywords are activated.  This should be tested.  Boundary
  ;; conditions can be wierd.
  if within then begin
      if numgti EQ 0 then $
        newgti = [ tstart, max(tt) ] $
      else $
        newgti = [ newgti, tstart, max(tt) ]
      numgti   = numgti + 1
      within   = 0
  endif
  
  ;; Empty GTI is a special case.  Here I have opted to *not* take the
  ;; lame HEASARC route.  Instead, I define a single element GTI array
  ;; to mean no good times.
  if numgti EQ 0 then begin
      EMPTY_GTI:
      count = 0L
      return, 0L
  endif

  ;; Reconstitute the new GTI array
  count    = numgti
  newgti   = reform(newgti, 2, numgti, /overwrite)

  ;; Now do some clean-up, removing zero-time GTIs and GTIs that
  ;; adjoin.  This really should not be necessary, since by definition
  ;; the above technique should not create buggy GTIs.
  if numgti GT 0 then begin
      ocount = count
      newgti = gtitrim(newgti, maxgap=ttol, count=count)
  endif

  return, newgti
end