source: trunk/SRC/Matrix/cmapply.pro @ 163

;+
; @file_comments
; Utility function, adapted from CMPRODUCT
;
; @param X
;
;
; @version $Id$
;
; @todo seb
;-
;
function cmapply_product, x
;
  compile_opt idl2, strictarrsubs
;
  sz = size(x)
  n = sz[1]

  while n GT 1 do begin
      if (n mod 2) EQ 1 then x[0,*] = x[0,*] * x[n-1,*]
      n2 = floor(n/2)
      x = x[0:n2-1,*] * x[n2:*,*]
      n = n2
  endwhile
  return, reform(x[0,*], /overwrite)
end

;+
; @file_comments
; cmapply_redim : Utility function, used to collect collaped dimensions
;
; @param newarr
;
;
; @param dimapply
;
;
; @param dimkeep
;
;
; @param nkeep
;
;
; @param totcol
;
;
; @param totkeep
;
;
; @todo seb
;-
;
pro cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
;
  compile_opt idl2, strictarrsubs
;
  sz = size(newarr)
  ;; First task: rearrange dimensions so that the dimensions
  ;; that are "kept" (ie, uncollapsed) are at the back
  dimkeep = where(histogram(dimapply,min=1,max=sz[0]) ne 1, nkeep)
  if nkeep EQ 0 then return

  newarr = transpose(temporary(newarr), [dimapply-1, dimkeep])
  ;; totcol is the total number of collapsed elements
  totcol = sz[dimapply[0]]
  for i = 1, n_elements(dimapply)-1 do totcol = totcol * sz[dimapply[i]]
  totkeep = sz[dimkeep[0]+1]
  for i = 1, n_elements(dimkeep)-1 do totkeep = totkeep * sz[dimkeep[i]+1]

  ;; this new array has two dimensions:
  ;;   * the first, all elements that will be collapsed
  ;;   * the second, all dimensions that will be preserved
  ;; (the ordering is so that all elements to be collapsed are
  ;;  adjacent in memory)
  newarr = reform(newarr, [totcol, totkeep], /overwrite)
end

;; Main function
;+
;
; @file_comments
; Applies a function to specified dimensions of an array
;
; Description:
;
; CMAPPLY will apply one of a few select functions to specified
; dimensions of an array.  Unlike some IDL functions, you *do* have
; a choice of which dimensions that are to be "collapsed" by this
; function.  Iterative loops are avoided where possible, for
; performance reasons.
;
;   The possible functions are:             (and number of loop iterations:)
;     +     - Performs a sum (as in TOTAL)       number of collapsed dimensions
;     AND   - Finds LOGICAL "AND" (not bitwise)  same
;     OR    - Finds LOGICAL "OR"  (not bitwise)  same
;     *     - Performs a product                 LOG_2[no. of collapsed elts.]
;
;     MIN   - Finds the minimum value            smaller of no. of collapsed
;     MAX   - Finds the maximum value            or output elements
;
;     USER  - Applies user-defined function      no. of output elements
;
;
;   It is possible to perform user-defined operations arrays using
;   CMAPPLY.  The OP parameter is set to 'USER:FUNCTNAME', where
;   FUNCTNAME is the name of a user-defined function.  The user
;   defined function should be defined such that it accepts a single
;   parameter, a vector, and returns a single scalar value.  Here is a
;   prototype for the function definition:
;
;      FUNCTION FUNCTNAME, x, KEYWORD1=key1, ...
;         scalar = ... function of x or keywords ...
;         RETURN, scalar
;      END
;
;   The function may accept keywords.  Keyword values are passed in to
;   CMAPPLY through the FUNCTARGS keywords parameter, and passed to
;   the user function via the _EXTRA mechanism.  Thus, while the
;   definition of the user function is highly constrained in the
;   number of positional parameters, there is absolute freedom in
;   passing keyword parameters.
;
;   It's worth noting however, that the implementation of user-defined
;   functions is not particularly optimized for speed.  Users are
;   encouraged to implement their own array if the number of output
;   elements is large.
;
; @categories 
; Array
;
; @param OP {in}{required}{type=string}
; The operation to perform, as a string.  May be upper or lower case.
;
; If a user-defined operation is to be passed, then OP is of
; the form, 'USER:FUNCTNAME', where FUNCTNAME is the name of
; the user-defined function.
;
; @param ARRAY {in}{required}{type=array} 
; An array of values to be operated on.
; Must not be of type STRING (7) or STRUCTURE (8).
;
; @param dimapply {in}{optional}{default=1 (ie, first dimension)}{type=array}
; An array of dimensions that are to be "collapsed", where
; the the first dimension starts with 1 (ie, same convention
; as IDL function TOTAL).  Whereas TOTAL only allows one
; dimension to be added, you can specify multiple dimensions
; to CMAPPLY.  Order does not matter, since all operations
; are associative and transitive.  NOTE: the dimensions refer
; to the *input* array, not the output array.  IDL allows a
; maximum of 8 dimensions.
;
; @keyword DOUBLE {default=not set}
; Set this if you wish the internal computations to be done
; in double precision if necessary.  If ARRAY is double
; precision (real or complex) then DOUBLE=1 is implied.
;
; @keyword TYPE {default=same as input type}
; Set this to the IDL code of the desired output type (refer
; to documentation of SIZE()).  Internal results will be
; rounded to the nearest integer if the output type is an
; integer type.
;
; @keyword FUNCTARGS
; If OP is 'USER:...', then the contents of this keyword
; are passed to the user function using the _EXTRA
; mechanism.  This way you can pass additional data to
; your user-supplied function, via keywords, without
; using common blocks.
; DEFAULT: undefined (i.e., no keywords passed by _EXTRA)
;
; @returns
; An array of the required TYPE, whose elements are the result of
; the requested operation.  Depending on the operation and number of
; elements in the input array, the result may be vulnerable to
; overflow or underflow.
;
; @examples
;
;   First example:  Shows how CMAPPLY can be used to total the second dimension of the
;   array called IN.  This is equivalent to OUT = TOTAL(IN, 2)
;
;   IDL> IN  = INDGEN(5,5)
;   IDL> OUT = CMAPPLY('+', IN, [2])
;   IDL> HELP, OUT
;   OUT             INT       = Array[5]
;
;   Second example:  Input is assumed to be an 5x100 array of 1's and
;   0's indicating the status of 5 detectors at 100 points in time.
;   The desired output is an array of 100 values, indicating whether
;   all 5 detectors are on (=1) at one time.  Use the logical AND
;   operation.
;
;   IDL> IN = detector_status             ; 5x100 array
;   IDL> OUT = CMAPPLY('AND', IN, [1])    ; collapses 1st dimension
;   IDL> HELP, OUT
;   OUT             BYTE      = Array[100]
;
;   (note that MIN could also have been used in this particular case,
;   although there would have been more loop iterations).
;
;   Third example:  Shows sum over first and third dimensions in an
;   array with dimensions 4x4x4:
;
;   IDL> IN = INDGEN(4,4,4)
;   IDL> OUT = CMAPPLY('+', IN, [1,3])
;   IDL> PRINT, OUT
;        408     472     536     600
;
;   Fourth example:  A user-function (MEDIAN) is used:
;
;   IDL> IN = RANDOMN(SEED,10,10,5)
;   IDL> OUT = CMAPPLY('USER:MEDIAN', IN, 3)
;   IDL> HELP, OUT
;   OUT             FLOAT     = Array[10, 10]
;
;   (OUT[i,j] is the median value of IN[i,j,*])
;
; @history Mar 1998, Written, CM
;   Changed usage message to not bomb, 24 Mar 2000, CM
;   Significant rewrite for *, MIN and MAX (inspired by Todd Clements
;     <Todd_Clements\@alumni.hmc.edu>); FOR loop indices are now type
;     LONG; copying terms are liberalized, CM, 22, Aug 2000
;   More efficient MAX/MIN (inspired by Alex Schuster), CM, 25 Jan
;     2002
;   Make new MAX/MIN actually work with 3d arrays, CM, 08 Feb 2002
;   Add user-defined functions, ON_ERROR, CM, 09 Feb 2002
;   Correct bug in MAX/MIN initialization of RESULT, CM, 05 Dec 2002
;
;  Author: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;  craigm\@lheamail.gsfc.nasa.gov
;
; @version $Id$
;
;-
function cmapply, op, array, dimapply, double=dbl, type=type, $
                  functargs=functargs, nocatch=nocatch
;
  compile_opt idl2, strictarrsubs
;

  if n_params() LT 2 then begin
      message, "USAGE: XX = CMAPPLY('OP',ARRAY,2)", /info
      message, '       where OP is +, *, AND, OR, MIN, MAX', /info
      return, -1L
  endif
  if NOT keyword_set(nocatch) then $
    on_error, 2 $
  else $
    on_error, 0

  ;; Parameter checking
  ;; 1) the dimensions of the array
  sz = size(array)
  if sz[0] EQ 0 then $
    message, 'ERROR: ARRAY must be an array!'

  ;; 2) The type of the array
  if sz[sz[0]+1] EQ 0 OR sz[sz[0]+1] EQ 7 OR sz[sz[0]+1] EQ 8 then $
    message, 'ERROR: Cannot apply to UNDEFINED, STRING, or STRUCTURE'
  if n_elements(type) EQ 0 then type = sz[sz[0]+1]

  ;; 3) The type of the operation
  szop = size(op)
  if szop[szop[0]+1] NE 7 then $
    message, 'ERROR: operation OP was not a string'

  ;; 4) The dimensions to apply (default is to apply to first dim)
  if n_params() EQ 2 then dimapply = 1
  dimapply = [ dimapply ]
  dimapply = dimapply[sort(dimapply)]   ; Sort in ascending order
  napply = n_elements(dimapply)

  ;; 5) Use double precision if requested or if needed
  if n_elements(dbl) EQ 0 then begin
      dbl=0
      if type EQ 5 OR type EQ 9 then dbl=1
  endif

  newop = strupcase(op)
  newarr = array
  newarr = reform(newarr, sz[1:sz[0]], /overwrite)
  case 1 of

      ;; *** Addition
      (newop EQ '+'): begin
          for i = 0L, napply-1 do begin
              newarr = total(temporary(newarr), dimapply[i]-i, double=dbl)
          endfor
      end

      ;; *** Multiplication
      (newop EQ '*'): begin ;; Multiplication (by summation of logarithms)
          cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
          if nkeep EQ 0 then begin
              newarr = reform(newarr, n_elements(newarr), 1, /overwrite)
              return, (cmapply_product(newarr))[0]
          endif

          result = cmapply_product(newarr)
          result = reform(result, sz[dimkeep+1], /overwrite)
          return, result
      end

      ;; *** LOGICAL AND or OR
      ((newop EQ 'AND') OR (newop EQ 'OR')): begin
          newarr = temporary(newarr) NE 0
          totelt = 1L
          for i = 0L, napply-1 do begin
              newarr = total(temporary(newarr), dimapply[i]-i)
              totelt = totelt * sz[dimapply[i]]
          endfor
          if newop EQ 'AND' then return, (round(newarr) EQ totelt)
          if newop EQ 'OR'  then return, (round(newarr) NE 0)
      end

      ;; Operations requiring a little more attention over how to
      ;; iterate
      ((newop EQ 'MAX') OR (newop EQ 'MIN')): begin
          cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
          if nkeep EQ 0 then begin
              if newop EQ 'MAX' then return, max(newarr)
              if newop EQ 'MIN' then return, min(newarr)
          endif

          ;; Next task: create result array
          result = make_array(totkeep, type=type)

          ;; Now either iterate over the number of output elements, or
          ;; the number of collapsed elements, whichever is smaller.
          if totcol LT totkeep then begin
              ;; Iterate over the number of collapsed elements
              result[0] = reform(newarr[0,*],totkeep,/overwrite)
              case newop of
                  'MAX': for i = 1L, totcol-1 do $
                    result[0] = result > newarr[i,*]
                  'MIN': for i = 1L, totcol-1 do $
                    result[0] = result < newarr[i,*]
              endcase
          endif else begin
              ;; Iterate over the number of output elements
              case newop of
                  'MAX': for i = 0L, totkeep-1 do result[i] = max(newarr[*,i])
                  'MIN': for i = 0L, totkeep-1 do result[i] = min(newarr[*,i])
              endcase
          endelse

          result = reform(result, sz[dimkeep+1], /overwrite)
          return, result
      end

      ;; User function
      (strmid(newop,0,4) EQ 'USER'): begin
          functname = strmid(newop,5)
          if functname EQ '' then $
            message, 'ERROR: '+newop+' is not a valid operation'

          cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
          if nkeep EQ 0 then begin
              if n_elements(functargs) GT 0 then $
                return, call_function(functname, newarr, _EXTRA=functargs)
              return, call_function(functname, newarr)
          endif

          ;; Next task: create result array
          result = make_array(totkeep, type=type)

          ;; Iterate over the number of output elements
          if n_elements(functargs) GT 0 then begin
              for i = 0L, totkeep-1 do $
                result[i] = call_function(functname, newarr[*,i], _EXTRA=functargs)
          endif else begin
              for i = 0L, totkeep-1 do $
                result[i] = call_function(functname, newarr[*,i])
          endelse

          result = reform(result, sz[dimkeep+1], /overwrite)
          return, result
      end


  endcase

  newsz = size(newarr)
  if type EQ newsz[newsz[0]+1] then return, newarr

  ;; Cast the result into the desired type, if necessary
  castfns = ['UNDEF', 'BYTE', 'FIX', 'LONG', 'FLOAT', $
             'DOUBLE', 'COMPLEX', 'UNDEF', 'UNDEF', 'DCOMPLEX' ]
  if type GE 1 AND type LE 3 then $
    return, call_function(castfns[type], round(newarr)) $
  else $
    return, call_function(castfns[type], newarr)
end