Changeset 31 for trunk/ToBeReviewed/MATRICE/cmapply.pro

Timestamp:

05/02/06 15:54:11 (18 years ago)

Author:

pinsard

Message:

upgrade of MATRICE according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files

File:

• trunk/ToBeReviewed/MATRICE/cmapply.pro (copied) (copied from trunk/MATRICE/cmapply.pro) (9 diffs)

trunk/ToBeReviewed/MATRICE/cmapply.pro

@@ -25 +25 @@
 ;   The possible functions are:             (and number of loop iterations:)
 ;     +     - Performs a sum (as in TOTAL)       number of collapsed dimensions
-;     *     - Performs a product                 same
 ;     AND   - Finds LOGICAL "AND" (not bitwise)  same
 ;     OR    - Finds LOGICAL "OR"  (not bitwise)  same
-;
-;     MIN   - Finds the minimum value            number of output elements
-;     MAX   - Finds the maximum value            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.
+;
 ;
 ; INPUTS:
@@ -36 +64 @@
 ;   OP - 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.
 ;
 ;   ARRAY - An array of values to be operated on.  Must not be of type
@@ -65 +97 @@
 ;          DEFAULT: same is input type
 ;
+;   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)
+;
 ; RETURN VALUE:
 ;
@@ -103 +142 @@
 ;        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,*))
+;
 ; MODIFICATION HISTORY:
 ;   Mar 1998, Written, CM
+;   Changed usage message to not bomb, 24 Mar 2000, CM
+;   Signficant 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
+;
+;  $Id$
 ;
 ;-
-
-function cmapply, op, array, dimapply, double=dbl, type=type
+; Copyright (C) 1998, 2000, 2002, 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.
+;-
+
+;; Utility function, adapted from CMPRODUCT
+function cmapply_product, x
+  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
+
+;; Utility function, used to collect collaped dimensions
+pro cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
+  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
+function cmapply, op, array, dimapply, double=dbl, type=type, $
+                  functargs=functargs, nocatch=nocatch
 
   if n_params() LT 2 then begin
       message, "USAGE: XX = CMAPPLY('OP',ARRAY,2)", /info
-      message, '       where OP is +, *, AND, OR, MIN, MAX'
+      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
@@ -151 +259 @@
       ;; *** Addition
       (newop EQ '+'): begin
-          for i = 0, napply-1 do begin
+          for i = 0L, napply-1 do begin
               newarr = total(temporary(newarr), dimapply(i)-i, double=dbl)
           endfor
@@ -158 +266 @@
       ;; *** Multiplication
       (newop EQ '*'): begin ;; Multiplication (by summation of logarithms)
-          dummy = check_math(1, 1)  ;; Disable printing of messages
-          ;; The following step seems to be the slowest
-          newarr = alog(abs(dcomplex(temporary(newarr))))
-          for i = 0, napply-1 do begin
-              newarr = total(temporary(newarr), dimapply(i)-i, double=1)
-          endfor
-          newarr = exp(temporary(newarr))
-          dummy = check_math(0, 0)  ;; Enable printing of messages
+          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
 
@@ -172 +281 @@
           newarr = temporary(newarr) NE 0
           totelt = 1L
-          for i = 0, napply-1 do begin
+          for i = 0L, napply-1 do begin
               newarr = total(temporary(newarr), dimapply(i)-i)
               totelt = totelt * sz(dimapply(i))
@@ -180 +289 @@
       end
 
-      ;; *** Operations requiring element by element access ... ho hum
+      ;; Operations requiring a little more attention over how to
+      ;; iterate
       ((newop EQ 'MAX') OR (newop EQ 'MIN')): begin
-          ;; 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, count)
-          if count EQ 0 then begin
-              case newop of
-                  'MAX': return, max(newarr)
-                  'MIN': return, min(newarr)
-              endcase
+          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
-          newarr = transpose( temporary(newarr), [ dimapply-1, dimkeep ])
-          ;; totcol is the total number of collapsed elements
-          totcol = exp(total(alog(sz(dimapply))))
-          totkeep = exp(total(alog(sz(dimkeep+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, round([totcol, totkeep]), /overwrite)
           
           ;; Next task: create result array
-          result = make_array(dimension=sz(dimkeep+1), type=type)
-
-          ;; Finally, compute the result, element by element
-          case newop of 
-              'MAX': for i = 0, totkeep-1 do result(i) = max(newarr(*,i))
-              'MIN': for i = 0, totkeep-1 do result(i) = min(newarr(*,i))
-          endcase
+          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
@@ -226 +367 @@
     return, call_function(castfns(type), newarr)
 end
+