;+
;
; @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 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
; ); 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