- Timestamp:
- 05/02/06 15:54:11 (18 years ago)
- Location:
- trunk
- Files:
-
- 1 added
- 1 deleted
- 9 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/ToBeReviewed/MATRICE/cmapply.pro
r29 r31 25 25 ; The possible functions are: (and number of loop iterations:) 26 26 ; + - Performs a sum (as in TOTAL) number of collapsed dimensions 27 ; * - Performs a product same28 27 ; AND - Finds LOGICAL "AND" (not bitwise) same 29 28 ; OR - Finds LOGICAL "OR" (not bitwise) same 30 ; 31 ; MIN - Finds the minimum value number of output elements 32 ; MAX - Finds the maximum value same 29 ; * - Performs a product LOG_2[no. of collapsed elts.] 30 ; 31 ; MIN - Finds the minimum value smaller of no. of collapsed 32 ; MAX - Finds the maximum value or output elements 33 ; 34 ; USER - Applies user-defined function no. of output elements 35 ; 36 ; 37 ; It is possible to perform user-defined operations arrays using 38 ; CMAPPLY. The OP parameter is set to 'USER:FUNCTNAME', where 39 ; FUNCTNAME is the name of a user-defined function. The user 40 ; defined function should be defined such that it accepts a single 41 ; parameter, a vector, and returns a single scalar value. Here is a 42 ; prototype for the function definition: 43 ; 44 ; FUNCTION FUNCTNAME, x, KEYWORD1=key1, ... 45 ; scalar = ... function of x or keywords ... 46 ; RETURN, scalar 47 ; END 48 ; 49 ; The function may accept keywords. Keyword values are passed in to 50 ; CMAPPLY through the FUNCTARGS keywords parameter, and passed to 51 ; the user function via the _EXTRA mechanism. Thus, while the 52 ; definition of the user function is highly constrained in the 53 ; number of positional parameters, there is absolute freedom in 54 ; passing keyword parameters. 55 ; 56 ; It's worth noting however, that the implementation of user-defined 57 ; functions is not particularly optimized for speed. Users are 58 ; encouraged to implement their own array if the number of output 59 ; elements is large. 60 ; 33 61 ; 34 62 ; INPUTS: … … 36 64 ; OP - The operation to perform, as a string. May be upper or lower 37 65 ; case. 66 ; 67 ; If a user-defined operation is to be passed, then OP is of 68 ; the form, 'USER:FUNCTNAME', where FUNCTNAME is the name of 69 ; the user-defined function. 38 70 ; 39 71 ; ARRAY - An array of values to be operated on. Must not be of type … … 65 97 ; DEFAULT: same is input type 66 98 ; 99 ; FUNCTARGS - If OP is 'USER:...', then the contents of this keyword 100 ; are passed to the user function using the _EXTRA 101 ; mechanism. This way you can pass additional data to 102 ; your user-supplied function, via keywords, without 103 ; using common blocks. 104 ; DEFAULT: undefined (i.e., no keywords passed by _EXTRA) 105 ; 67 106 ; RETURN VALUE: 68 107 ; … … 103 142 ; 408 472 536 600 104 143 ; 144 ; Fourth example. A user-function (MEDIAN) is used: 145 ; 146 ; IDL> IN = RANDOMN(SEED,10,10,5) 147 ; IDL> OUT = CMAPPLY('USER:MEDIAN', IN, 3) 148 ; IDL> HELP, OUT 149 ; OUT FLOAT = Array[10, 10] 150 ; 151 ; (OUT(i,j) is the median value of IN(i,j,*)) 152 ; 105 153 ; MODIFICATION HISTORY: 106 154 ; Mar 1998, Written, CM 155 ; Changed usage message to not bomb, 24 Mar 2000, CM 156 ; Signficant rewrite for *, MIN and MAX (inspired by Todd Clements 157 ; <Todd_Clements@alumni.hmc.edu>); FOR loop indices are now type 158 ; LONG; copying terms are liberalized, CM, 22, Aug 2000 159 ; More efficient MAX/MIN (inspired by Alex Schuster), CM, 25 Jan 160 ; 2002 161 ; Make new MAX/MIN actually work with 3d arrays, CM, 08 Feb 2002 162 ; Add user-defined functions, ON_ERROR, CM, 09 Feb 2002 163 ; Correct bug in MAX/MIN initialization of RESULT, CM, 05 Dec 2002 164 ; 165 ; $Id$ 107 166 ; 108 167 ;- 109 110 function cmapply, op, array, dimapply, double=dbl, type=type 168 ; Copyright (C) 1998, 2000, 2002, Craig Markwardt 169 ; This software is provided as is without any warranty whatsoever. 170 ; Permission to use, copy, modify, and distribute modified or 171 ; unmodified copies is granted, provided this copyright and disclaimer 172 ; are included unchanged. 173 ;- 174 175 ;; Utility function, adapted from CMPRODUCT 176 function cmapply_product, x 177 sz = size(x) 178 n = sz(1) 179 180 while n GT 1 do begin 181 if (n mod 2) EQ 1 then x(0,*) = x(0,*) * x(n-1,*) 182 n2 = floor(n/2) 183 x = x(0:n2-1,*) * x(n2:*,*) 184 n = n2 185 endwhile 186 return, reform(x(0,*), /overwrite) 187 end 188 189 ;; Utility function, used to collect collaped dimensions 190 pro cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep 191 sz = size(newarr) 192 ;; First task: rearrange dimensions so that the dimensions 193 ;; that are "kept" (ie, uncollapsed) are at the back 194 dimkeep = where(histogram(dimapply,min=1,max=sz(0)) ne 1, nkeep) 195 if nkeep EQ 0 then return 196 197 newarr = transpose(temporary(newarr), [dimapply-1, dimkeep]) 198 ;; totcol is the total number of collapsed elements 199 totcol = sz(dimapply(0)) 200 for i = 1, n_elements(dimapply)-1 do totcol = totcol * sz(dimapply(i)) 201 totkeep = sz(dimkeep(0)+1) 202 for i = 1, n_elements(dimkeep)-1 do totkeep = totkeep * sz(dimkeep(i)+1) 203 204 ;; this new array has two dimensions: 205 ;; * the first, all elements that will be collapsed 206 ;; * the second, all dimensions that will be preserved 207 ;; (the ordering is so that all elements to be collapsed are 208 ;; adjacent in memory) 209 newarr = reform(newarr, [totcol, totkeep], /overwrite) 210 end 211 212 ;; Main function 213 function cmapply, op, array, dimapply, double=dbl, type=type, $ 214 functargs=functargs, nocatch=nocatch 111 215 112 216 if n_params() LT 2 then begin 113 217 message, "USAGE: XX = CMAPPLY('OP',ARRAY,2)", /info 114 message, ' where OP is +, *, AND, OR, MIN, MAX' 218 message, ' where OP is +, *, AND, OR, MIN, MAX', /info 115 219 return, -1L 116 220 endif 221 if NOT keyword_set(nocatch) then $ 222 on_error, 2 $ 223 else $ 224 on_error, 0 117 225 118 226 ;; Parameter checking … … 151 259 ;; *** Addition 152 260 (newop EQ '+'): begin 153 for i = 0 , napply-1 do begin261 for i = 0L, napply-1 do begin 154 262 newarr = total(temporary(newarr), dimapply(i)-i, double=dbl) 155 263 endfor … … 158 266 ;; *** Multiplication 159 267 (newop EQ '*'): begin ;; Multiplication (by summation of logarithms) 160 dummy = check_math(1, 1) ;; Disable printing of messages 161 ;; The following step seems to be the slowest 162 newarr = alog(abs(dcomplex(temporary(newarr)))) 163 for i = 0, napply-1 do begin 164 newarr = total(temporary(newarr), dimapply(i)-i, double=1) 165 endfor 166 newarr = exp(temporary(newarr)) 167 dummy = check_math(0, 0) ;; Enable printing of messages 268 cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep 269 if nkeep EQ 0 then begin 270 newarr = reform(newarr, n_elements(newarr), 1, /overwrite) 271 return, (cmapply_product(newarr))(0) 272 endif 273 274 result = cmapply_product(newarr) 275 result = reform(result, sz(dimkeep+1), /overwrite) 276 return, result 168 277 end 169 278 … … 172 281 newarr = temporary(newarr) NE 0 173 282 totelt = 1L 174 for i = 0 , napply-1 do begin283 for i = 0L, napply-1 do begin 175 284 newarr = total(temporary(newarr), dimapply(i)-i) 176 285 totelt = totelt * sz(dimapply(i)) … … 180 289 end 181 290 182 ;; *** Operations requiring element by element access ... ho hum 291 ;; Operations requiring a little more attention over how to 292 ;; iterate 183 293 ((newop EQ 'MAX') OR (newop EQ 'MIN')): begin 184 ;; First task: rearrange dimensions so that the dimensions 185 ;; that are "kept" (ie, uncollapsed) are at the back 186 dimkeep = where(histogram(dimapply,min=1,max=sz(0)) ne 1, count) 187 if count EQ 0 then begin 188 case newop of 189 'MAX': return, max(newarr) 190 'MIN': return, min(newarr) 191 endcase 294 cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep 295 if nkeep EQ 0 then begin 296 if newop EQ 'MAX' then return, max(newarr) 297 if newop EQ 'MIN' then return, min(newarr) 192 298 endif 193 newarr = transpose( temporary(newarr), [ dimapply-1, dimkeep ])194 ;; totcol is the total number of collapsed elements195 totcol = exp(total(alog(sz(dimapply))))196 totkeep = exp(total(alog(sz(dimkeep+1))))197 ;; this new array has two dimensions:198 ;; * the first, all elements that will be collapsed199 ;; * the second, all dimensions that will be preserved200 ;; (the ordering is so that all elements to be collapsed are201 ;; adjacent in memory)202 newarr = reform(newarr, round([totcol, totkeep]), /overwrite)203 299 204 300 ;; Next task: create result array 205 result = make_array(dimension=sz(dimkeep+1), type=type) 206 207 ;; Finally, compute the result, element by element 208 case newop of 209 'MAX': for i = 0, totkeep-1 do result(i) = max(newarr(*,i)) 210 'MIN': for i = 0, totkeep-1 do result(i) = min(newarr(*,i)) 211 endcase 301 result = make_array(totkeep, type=type) 302 303 ;; Now either iterate over the number of output elements, or 304 ;; the number of collapsed elements, whichever is smaller. 305 if totcol LT totkeep then begin 306 ;; Iterate over the number of collapsed elements 307 result(0) = reform(newarr(0,*),totkeep,/overwrite) 308 case newop of 309 'MAX': for i = 1L, totcol-1 do $ 310 result(0) = result > newarr(i,*) 311 'MIN': for i = 1L, totcol-1 do $ 312 result(0) = result < newarr(i,*) 313 endcase 314 endif else begin 315 ;; Iterate over the number of output elements 316 case newop of 317 'MAX': for i = 0L, totkeep-1 do result(i) = max(newarr(*,i)) 318 'MIN': for i = 0L, totkeep-1 do result(i) = min(newarr(*,i)) 319 endcase 320 endelse 321 322 result = reform(result, sz(dimkeep+1), /overwrite) 212 323 return, result 213 324 end 325 326 ;; User function 327 (strmid(newop,0,4) EQ 'USER'): begin 328 functname = strmid(newop,5) 329 if functname EQ '' then $ 330 message, 'ERROR: '+newop+' is not a valid operation' 331 332 cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep 333 if nkeep EQ 0 then begin 334 if n_elements(functargs) GT 0 then $ 335 return, call_function(functname, newarr, _EXTRA=functargs) 336 return, call_function(functname, newarr) 337 endif 338 339 ;; Next task: create result array 340 result = make_array(totkeep, type=type) 341 342 ;; Iterate over the number of output elements 343 if n_elements(functargs) GT 0 then begin 344 for i = 0L, totkeep-1 do $ 345 result(i) = call_function(functname, newarr(*,i), _EXTRA=functargs) 346 endif else begin 347 for i = 0L, totkeep-1 do $ 348 result(i) = call_function(functname, newarr(*,i)) 349 endelse 350 351 result = reform(result, sz(dimkeep+1), /overwrite) 352 return, result 353 end 354 214 355 215 356 endcase … … 226 367 return, call_function(castfns(type), newarr) 227 368 end 369 -
trunk/ToBeReviewed/MATRICE/colle.pro
r29 r31 133 133 ; pour cela on contient le vecteur permute qui donne la place que 134 134 ; doivent prendre les dimensions ds la matrice transposee 135 siz = (size(*ptrtab[0]))[0] 136 if siz LT direc then $ 137 *ptrtab[0] = reform(*ptrtab[0], [(size(*ptrtab[0]))[1:siz], replicate(1, direc-siz)], /over) 135 138 permute = indgen((size(*ptrtab[0]))[0]) 136 139 permute[0] = direc-1 137 140 permute[direc-1] = 0 138 if (size(*ptrtab[0]))[0] NE direc then $139 *ptrtab[0] = reform(*ptrtab[0], [(size(*ptrtab[0]))[1:direc-1], 1], /over)140 141 res = transpose(*ptrtab[0], permute) 141 142 if NOT keyword_set(sauve) then ptr_free, ptrtab[0] 142 143 FOR n = 1,nbretab-1 DO BEGIN ; on colle suivant la dimension 1 143 if (size(*ptrtab[n]))[0] NEdirec then $144 *ptrtab[n] = reform(*ptrtab[n], [(size(*ptrtab[n]))[1: direc-1], 1])144 if (size(*ptrtab[n]))[0] LT direc then $ 145 *ptrtab[n] = reform(*ptrtab[n], [(size(*ptrtab[n]))[1:siz], replicate(1, direc-siz)]) 145 146 res = [temporary(res), transpose(*ptrtab[n], permute)] 146 147 if NOT keyword_set(sauve) then ptr_free, ptrtab[n] -
trunk/ToBeReviewed/MATRICE/inter.pro
r29 r31 11 11 ; CALLING SEQUENCE:res=inter(a,b) 12 12 ; 13 ; INPUTS:a et b:arrays of positive integers, which need 14 ; not be sorted. Duplicate elements are ignored, as they have no15 ; effect on theresult13 ; INPUTS:a et b:arrays of positive integers, which need not to be 14 ; sorted. Duplicate elements are ignored, as they have noeffect on the 15 ; result 16 16 ; 17 17 ; KEYWORD PARAMETERS:
Note: See TracChangeset
for help on using the changeset viewer.