1Primitive(3) User Contributed Perl Documentation Primitive(3)
2
6 PDL::Primitive - primitive operations for pdl
7
9 This module provides some primitive and useful functions defined using
10 PDL::PP and able to use the new indexing tricks.
11 See PDL::Indexing for how to use indices creatively. For explanation
13 of the signature format, see PDL::PP.
14
16 # Pulls in PDL::Primitive, among other modules.
17 use PDL;
18 # Only pull in PDL::Primitive:
20 use PDL::Primitive;
21
23 inner
24 Signature: (a(n); b(n); [o]c())
25 Inner product over one dimension
27 c = sum_i a_i * b_i
29 If "a() * b()" contains only bad data, "c()" is set bad. Otherwise
31 "c()" will have its bad flag cleared, as it will not contain any bad
32 values.
33 outer
35 Signature: (a(n); b(m); [o]c(n,m))
36 outer product over one dimension
38 Naturally, it is possible to achieve the effects of outer product
40 simply by broadcasting over the ""*"" operator but this function is
41 provided for convenience.
42 outer processes bad values. It will set the bad-value flag of all
44 output ndarrays if the flag is set for any of the input ndarrays.
45 x
47 Signature: (a(i,z), b(x,i),[o]c(x,z))
48 Matrix multiplication
50 PDL overloads the "x" operator (normally the repeat operator) for
52 matrix multiplication. The number of columns (size of the 0 dimension)
53 in the left-hand argument must normally equal the number of rows (size
54 of the 1 dimension) in the right-hand argument.
55 Row vectors are represented as (N x 1) two-dimensional PDLs, or you may
57 be sloppy and use a one-dimensional PDL. Column vectors are
58 represented as (1 x N) two-dimensional PDLs.
59 Broadcasting occurs in the usual way, but as both the 0 and 1 dimension
61 (if present) are included in the operation, you must be sure that you
62 don't try to broadcast over either of those dims.
63 Of note, due to how Perl v5.14.0 and above implement operator
65 overloading of the "x" operator, the use of parentheses for the left
66 operand creates a list context, that is
67 pdl> ( $x * $y ) x $z
69 ERROR: Argument "..." isn't numeric in repeat (x) ...
70 treats $z as a numeric count for the list repeat operation and does not
72 call the scalar form of the overloaded operator. To use the operator in
73 this case, use a scalar context:
74 pdl> scalar( $x * $y ) x $z
76 or by calling "matmult" directly:
78 pdl> ( $x * $y )->matmult( $z )
80 EXAMPLES
82 Here are some simple ways to define vectors and matrices:
84 pdl> $r = pdl(1,2); # A row vector
86 pdl> $c = pdl([[3],[4]]); # A column vector
87 pdl> $c = pdl(3,4)->(*1); # A column vector, using NiceSlice
88 pdl> $m = pdl([[1,2],[3,4]]); # A 2x2 matrix
89 Now that we have a few objects prepared, here is how to matrix-multiply
91 them:
92 pdl> print $r x $m # row x matrix = row
94 [
95 [ 7 10]
96 ]
97 pdl> print $m x $r # matrix x row = ERROR
99 PDL: Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1
100 pdl> print $m x $c # matrix x column = column
102 [
103 [ 5]
104 [11]
105 ]
106 pdl> print $m x 2 # Trivial case: scalar mult.
108 [
109 [2 4]
110 [6 8]
111 ]
112 pdl> print $r x $c # row x column = scalar
114 [
115 [11]
116 ]
117 pdl> print $c x $r # column x row = matrix
119 [
120 [3 6]
121 [4 8]
122 ]
123 INTERNALS
125 The mechanics of the multiplication are carried out by the "matmult"
127 method.
128 matmult
130 Signature: (a(t,h); b(w,t); [o]c(w,h))
131 Matrix multiplication
133 Notionally, matrix multiplication $x x $y is equivalent to the
135 broadcasting expression
136 $x->dummy(1)->inner($y->xchg(0,1)->dummy(2),$c);
138 but for large matrices that breaks CPU cache and is slow. Instead,
140 matmult calculates its result in 32x32x32 tiles, to keep the memory
141 footprint within cache as long as possible on most modern CPUs.
142 For usage, see "x", a description of the overloaded 'x' operator
144 matmult ignores the bad-value flag of the input ndarrays. It will set
146 the bad-value flag of all output ndarrays if the flag is set for any of
147 the input ndarrays.
148 innerwt
150 Signature: (a(n); b(n); c(n); [o]d())
151 Weighted (i.e. triple) inner product
153 d = sum_i a(i) b(i) c(i)
155 innerwt processes bad values. It will set the bad-value flag of all
157 output ndarrays if the flag is set for any of the input ndarrays.
158 inner2
160 Signature: (a(n); b(n,m); c(m); [o]d())
161 Inner product of two vectors and a matrix
163 d = sum_ij a(i) b(i,j) c(j)
165 Note that you should probably not broadcast over "a" and "c" since that
167 would be very wasteful. Instead, you should use a temporary for "b*c".
168 inner2 processes bad values. It will set the bad-value flag of all
170 output ndarrays if the flag is set for any of the input ndarrays.
171 inner2d
173 Signature: (a(n,m); b(n,m); [o]c())
174 Inner product over 2 dimensions.
176 Equivalent to
178 $c = inner($x->clump(2), $y->clump(2))
180 inner2d processes bad values. It will set the bad-value flag of all
182 output ndarrays if the flag is set for any of the input ndarrays.
183 inner2t
185 Signature: (a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k)))
186 Efficient Triple matrix product "a*b*c"
188 Efficiency comes from by using the temporary "tmp". This operation only
190 scales as "N**3" whereas broadcasting using "inner2" would scale as
191 "N**4".
192 The reason for having this routine is that you do not need to have the
194 same broadcast-dimensions for "tmp" as for the other arguments, which
195 in case of large numbers of matrices makes this much more memory-
196 efficient.
197 It is hoped that things like this could be taken care of as a kind of
199 closures at some point.
200 inner2t processes bad values. It will set the bad-value flag of all
202 output ndarrays if the flag is set for any of the input ndarrays.
203 crossp
205 Signature: (a(tri=3); b(tri); [o] c(tri))
206 Cross product of two 3D vectors
208 After
210 $c = crossp $x, $y
212 the inner product "$c*$x" and "$c*$y" will be zero, i.e. $c is
214 orthogonal to $x and $y
215 crossp does not process bad values. It will set the bad-value flag of
217 all output ndarrays if the flag is set for any of the input ndarrays.
218 norm
220 Signature: (vec(n); [o] norm(n))
221 Normalises a vector to unit Euclidean length
223 norm processes bad values. It will set the bad-value flag of all
225 output ndarrays if the flag is set for any of the input ndarrays.
226 indadd
228 Signature: (input(n); indx ind(n); [io] sum(m))
229 Broadcasting index add: add "input" to the "ind" element of "sum", i.e:
231 sum(ind) += input
233 Simple example:
235 $x = 2;
237 $ind = 3;
238 $sum = zeroes(10);
239 indadd($x,$ind, $sum);
240 print $sum
241 #Result: ( 2 added to element 3 of $sum)
242 # [0 0 0 2 0 0 0 0 0 0]
243 Broadcasting example:
245 $x = pdl( 1,2,3);
247 $ind = pdl( 1,4,6);
248 $sum = zeroes(10);
249 indadd($x,$ind, $sum);
250 print $sum."\n";
251 #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
252 # [0 1 0 0 2 0 3 0 0 0]
253 The routine barfs on bad indices, and bad inputs set target outputs
255 bad.
256 conv1d
258 Signature: (a(m); kern(p); [o]b(m); int reflect)
259 1D convolution along first dimension
261 The m-th element of the discrete convolution of an input ndarray $a of
263 size $M, and a kernel ndarray $kern of size $P, is calculated as
264 n = ($P-1)/2
266 ====
267 \
268 ($a conv1d $kern)[m] = > $a_ext[m - n] * $kern[n]
269 /
270 ====
271 n = -($P-1)/2
272 where $a_ext is either the periodic (or reflected) extension of $a so
274 it is equal to $a on " 0..$M-1 " and equal to the corresponding
275 periodic/reflected image of $a outside that range.
276 $con = conv1d sequence(10), pdl(-1,0,1);
278 $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
280 By default, periodic boundary conditions are assumed (i.e. wrap
282 around). Alternatively, you can request reflective boundary conditions
283 using the "Boundary" option:
284 {Boundary => 'reflect'} # case in 'reflect' doesn't matter
286 The convolution is performed along the first dimension. To apply it
288 across another dimension use the slicing routines, e.g.
289 $y = $x->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
291 This function is useful for broadcasted filtering of 1D signals.
293 Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv
295 WARNING: "conv1d" processes bad values in its inputs as the numeric
297 value of "$pdl->badvalue" so it is not recommended for processing pdls
298 with bad values in them unless special care is taken.
299 conv1d ignores the bad-value flag of the input ndarrays. It will set
301 the bad-value flag of all output ndarrays if the flag is set for any of
302 the input ndarrays.
303 in
305 Signature: (a(); b(n); [o] c())
306 test if a is in the set of values b
308 $goodmsk = $labels->in($goodlabels);
310 print pdl(3,1,4,6,2)->in(pdl(2,3,3));
311 [1 0 0 0 1]
312 "in" is akin to the is an element of of set theory. In principle, PDL
314 broadcasting could be used to achieve its functionality by using a
315 construct like
316 $msk = ($labels->dummy(0) == $goodlabels)->orover;
318 However, "in" doesn't create a (potentially large) intermediate and is
320 generally faster.
321 in does not process bad values. It will set the bad-value flag of all
323 output ndarrays if the flag is set for any of the input ndarrays.
324 uniq
326 return all unique elements of an ndarray
327 The unique elements are returned in ascending order.
329 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniq
331 [-1 0 2 4 6] # 0 is returned 2nd (sorted order)
332 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniq
334 [-1 2 4 6 nan] # NaN value is returned at end
335 Note: The returned pdl is 1D; any structure of the input ndarray is
337 lost. "NaN" values are never compare equal to any other values, even
338 themselves. As a result, they are always unique. "uniq" returns the
339 NaN values at the end of the result ndarray. This follows the Matlab
340 usage.
341 See "uniqind" if you need the indices of the unique elements rather
343 than the values.
344 Bad values are not considered unique by uniq and are ignored.
346 $x=sequence(10);
348 $x=$x->setbadif($x%3);
349 print $x->uniq;
350 [0 3 6 9]
351 uniqind
353 Return the indices of all unique elements of an ndarray The order is in
354 the order of the values to be consistent with uniq. "NaN" values never
355 compare equal with any other value and so are always unique. This
356 follows the Matlab usage.
357 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniqind
359 [5 4 1 3 6] # the 0 at index 4 is returned 2nd, but...
360 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniqind
362 [5 1 3 6 4] # ...the NaN at index 4 is returned at end
363 Note: The returned pdl is 1D; any structure of the input ndarray is
365 lost.
366 See "uniq" if you want the unique values instead of the indices.
368 Bad values are not considered unique by uniqind and are ignored.
370 uniqvec
372 Return all unique vectors out of a collection
373 NOTE: If any vectors in the input ndarray have NaN values
375 they are returned at the end of the non-NaN ones. This is
376 because, by definition, NaN values never compare equal with
377 any other value.
378 NOTE: The current implementation does not sort the vectors
380 containing NaN values.
381 The unique vectors are returned in lexicographically sorted ascending
383 order. The 0th dimension of the input PDL is treated as a dimensional
384 index within each vector, and the 1st and any higher dimensions are
385 taken to run across vectors. The return value is always 2D; any
386 structure of the input PDL (beyond using the 0th dimension for vector
387 index) is lost.
388 See also "uniq" for a unique list of scalars; and qsortvec for sorting
390 a list of vectors lexicographcally.
391 If a vector contains all bad values, it is ignored as in "uniq". If
393 some of the values are good, it is treated as a normal vector. For
394 example, [1 2 BAD] and [BAD 2 3] could be returned, but [BAD BAD BAD]
395 could not. Vectors containing BAD values will be returned after any
396 non-NaN and non-BAD containing vectors, followed by the NaN vectors.
397 hclip
399 Signature: (a(); b(); [o] c())
400 clip (threshold) $a by $b ($b is upper bound)
402 hclip processes bad values. It will set the bad-value flag of all
404 output ndarrays if the flag is set for any of the input ndarrays.
405 lclip
407 Signature: (a(); b(); [o] c())
408 clip (threshold) $a by $b ($b is lower bound)
410 lclip processes bad values. It will set the bad-value flag of all
412 output ndarrays if the flag is set for any of the input ndarrays.
413 clip
415 Clip (threshold) an ndarray by (optional) upper or lower bounds.
416 $y = $x->clip(0,3);
418 $c = $x->clip(undef, $x);
419 clip handles bad values since it is just a wrapper around "hclip" and
421 "lclip".
422 clip
424 Signature: (a(); l(); h(); [o] c())
425 info not available
427 clip processes bad values. It will set the bad-value flag of all
429 output ndarrays if the flag is set for any of the input ndarrays.
430 wtstat
432 Signature: (a(n); wt(n); avg(); [o]b(); int deg)
433 Weighted statistical moment of given degree
435 This calculates a weighted statistic over the vector "a". The formula
437 is
438 b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
440 Bad values are ignored in any calculation; $b will only have its bad
442 flag set if the output contains any bad data.
443 statsover
445 Signature: (a(n); w(n); float+ [o]avg(); float+ [o]prms(); int+ [o]median(); int+ [o]min(); int+ [o]max(); float+ [o]adev(); float+ [o]rms())
446 Calculate useful statistics over a dimension of an ndarray
448 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($ndarray, $weights);
450 This utility function calculates various useful quantities of an
452 ndarray. These are:
453 • the mean:
455 MEAN = sum (x)/ N
457 with "N" being the number of elements in x
459 • the population RMS deviation from the mean:
461 PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
463 The population deviation is the best-estimate of the deviation of
465 the population from which a sample is drawn.
466 • the median
468 The median is the 50th percentile data value. Median is found by
470 medover, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
471 • the minimum
473 • the maximum
475 • the average absolute deviation:
477 AADEV = sum( abs(x-mean(x)) )/N
479 • RMS deviation from the mean:
481 RMS = sqrt(sum( (x-mean(x))^2 )/N)
483 (also known as the root-mean-square deviation, or the square root of
485 the variance)
486 This operator is a projection operator so the calculation will take
488 place over the final dimension. Thus if the input is N-dimensional each
489 returned value will be N-1 dimensional, to calculate the statistics for
490 the entire ndarray either use "clump(-1)" directly on the ndarray or
491 call "stats".
492 Bad values are simply ignored in the calculation, effectively reducing
494 the sample size. If all data are bad then the output data are marked
495 bad.
496 stats
498 Calculates useful statistics on an ndarray
499 ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($ndarray,[$weights]);
501 This utility calculates all the most useful quantities in one call. It
503 works the same way as "statsover", except that the quantities are
504 calculated considering the entire input PDL as a single sample, rather
505 than as a collection of rows. See "statsover" for definitions of the
506 returned quantities.
507 Bad values are handled; if all input values are bad, then all of the
509 output values are flagged bad.
510 histogram
512 Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
513 Calculates a histogram for given stepsize and minimum.
515 $h = histogram($data, $step, $min, $numbins);
517 $hist = zeroes $numbins; # Put histogram in existing ndarray.
518 histogram($data, $hist, $step, $min, $numbins);
519 The histogram will contain $numbins bins starting from $min, each $step
521 wide. The value in each bin is the number of values in $data that lie
522 within the bin limits.
523 Data below the lower limit is put in the first bin, and data above the
525 upper limit is put in the last bin.
526 The output is reset in a different broadcastloop so that you can take a
528 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
529 For a higher-level interface, see hist.
531 pdl> p histogram(pdl(1,1,2),1,0,3)
533 [0 2 1]
534 histogram processes bad values. It will set the bad-value flag of all
536 output ndarrays if the flag is set for any of the input ndarrays.
537 whistogram
539 Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
540 Calculates a histogram from weighted data for given stepsize and
542 minimum.
543 $h = whistogram($data, $weights, $step, $min, $numbins);
545 $hist = zeroes $numbins; # Put histogram in existing ndarray.
546 whistogram($data, $weights, $hist, $step, $min, $numbins);
547 The histogram will contain $numbins bins starting from $min, each $step
549 wide. The value in each bin is the sum of the values in $weights that
550 correspond to values in $data that lie within the bin limits.
551 Data below the lower limit is put in the first bin, and data above the
553 upper limit is put in the last bin.
554 The output is reset in a different broadcastloop so that you can take a
556 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
557 pdl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
559 [0 0.2 0.5 0]
560 whistogram processes bad values. It will set the bad-value flag of all
562 output ndarrays if the flag is set for any of the input ndarrays.
563 histogram2d
565 Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
566 double stepb; double minb; int mbsize => mb;)
567 Calculates a 2d histogram.
569 $h = histogram2d($datax, $datay, $stepx, $minx,
571 $nbinx, $stepy, $miny, $nbiny);
572 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing ndarray.
573 histogram2d($datax, $datay, $hist, $stepx, $minx,
574 $nbinx, $stepy, $miny, $nbiny);
575 The histogram will contain $nbinx x $nbiny bins, with the lower limits
577 of the first one at "($minx, $miny)", and with bin size "($stepx,
578 $stepy)". The value in each bin is the number of values in $datax and
579 $datay that lie within the bin limits.
580 Data below the lower limit is put in the first bin, and data above the
582 upper limit is put in the last bin.
583 pdl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
585 [
586 [0 0 0]
587 [0 2 2]
588 [0 1 0]
589 ]
590 histogram2d processes bad values. It will set the bad-value flag of
592 all output ndarrays if the flag is set for any of the input ndarrays.
593 whistogram2d
595 Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
596 double stepb; double minb; int mbsize => mb;)
597 Calculates a 2d histogram from weighted data.
599 $h = whistogram2d($datax, $datay, $weights,
601 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
602 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing ndarray.
603 whistogram2d($datax, $datay, $weights, $hist,
604 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
605 The histogram will contain $nbinx x $nbiny bins, with the lower limits
607 of the first one at "($minx, $miny)", and with bin size "($stepx,
608 $stepy)". The value in each bin is the sum of the values in $weights
609 that correspond to values in $datax and $datay that lie within the bin
610 limits.
611 Data below the lower limit is put in the first bin, and data above the
613 upper limit is put in the last bin.
614 pdl> p whistogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),pdl(0.1,0.2,0.3,0.4,0.5),1,0,3,1,0,3)
616 [
617 [ 0 0 0]
618 [ 0 0.5 0.9]
619 [ 0 0.1 0]
620 ]
621 whistogram2d processes bad values. It will set the bad-value flag of
623 all output ndarrays if the flag is set for any of the input ndarrays.
624 fibonacci
626 Signature: (i(n); indx [o]x(n))
627 Constructor - a vector with Fibonacci's sequence
629 fibonacci does not process bad values. It will set the bad-value flag
631 of all output ndarrays if the flag is set for any of the input
632 ndarrays.
633 append
635 Signature: (a(n); b(m); [o] c(mn))
636 append two ndarrays by concatenating along their first dimensions
638 $x = ones(2,4,7);
640 $y = sequence 5;
641 $c = $x->append($y); # size of $c is now (7,4,7) (a jumbo-ndarray ;)
642 "append" appends two ndarrays along their first dimensions. The rest of
644 the dimensions must be compatible in the broadcasting sense. The
645 resulting size of the first dimension is the sum of the sizes of the
646 first dimensions of the two argument ndarrays - i.e. "n + m".
647 Similar functions include "glue" (below), which can append more than
649 two ndarrays along an arbitrary dimension, and cat, which can append
650 more than two ndarrays that all have the same sized dimensions.
651 append does not process bad values. It will set the bad-value flag of
653 all output ndarrays if the flag is set for any of the input ndarrays.
654 glue
656 $c = $x->glue(<dim>,$y,...)
657 Glue two or more PDLs together along an arbitrary dimension (N-D
659 "append").
660 Sticks $x, $y, and all following arguments together along the specified
662 dimension. All other dimensions must be compatible in the broadcasting
663 sense.
664 Glue is permissive, in the sense that every PDL is treated as having an
666 infinite number of trivial dimensions of order 1 -- so "$x->glue(3,$y)"
667 works, even if $x and $y are only one dimensional.
668 If one of the PDLs has no elements, it is ignored. Likewise, if one of
670 them is actually the undefined value, it is treated as if it had no
671 elements.
672 If the first parameter is a defined perl scalar rather than a pdl, then
674 it is taken as a dimension along which to glue everything else, so you
675 can say "$cube = PDL::glue(3,@image_list);" if you like.
676 "glue" is implemented in pdl, using a combination of xchg and "append".
678 It should probably be updated (one day) to a pure PP function.
679 Similar functions include "append" (above), which appends only two
681 ndarrays along their first dimension, and cat, which can append more
682 than two ndarrays that all have the same sized dimensions.
683 cmpvec
685 Signature: (a(n); b(n); sbyte [o]c())
686 Compare two vectors lexicographically.
688 Returns -1 if a is less, 1 if greater, 0 if equal.
690 The output is bad if any input values up to the point of inequality are
692 bad - any after are ignored.
693 eqvec
695 Signature: (a(n); b(n); sbyte [o]c())
696 Compare two vectors, returning 1 if equal, 0 if not equal.
698 The output is bad if any input values are bad.
700 enumvec
702 Signature: (v(M,N); indx [o]k(N))
703 Enumerate a list of vectors with locally unique keys.
705 Given a sorted list of vectors $v, generate a vector $k containing
707 locally unique keys for the elements of $v (where an "element" is a
708 vector of length $M ocurring in $v).
709 Note that the keys returned in $k are only unique over a run of a
711 single vector in $v, so that each unique vector in $v has at least one
712 0 (zero) index in $k associated with it. If you need global keys, see
713 enumvecg().
714 Contributed by Bryan Jurish <moocow@cpan.org>.
716 enumvec does not process bad values. It will set the bad-value flag of
718 all output ndarrays if the flag is set for any of the input ndarrays.
719 enumvecg
721 Signature: (v(M,N); indx [o]k(N))
722 Enumerate a list of vectors with globally unique keys.
724 Given a sorted list of vectors $v, generate a vector $k containing
726 globally unique keys for the elements of $v (where an "element" is a
727 vector of length $M ocurring in $v). Basically does the same thing as:
728 $k = $v->vsearchvec($v->uniqvec);
730 ... but somewhat more efficiently.
732 Contributed by Bryan Jurish <moocow@cpan.org>.
734 enumvecg does not process bad values. It will set the bad-value flag
736 of all output ndarrays if the flag is set for any of the input
737 ndarrays.
738 vsearchvec
740 Signature: (find(M); which(M,N); indx [o]found())
741 Routine for searching N-dimensional values - akin to vsearch() for
743 vectors.
744 $found = vsearchvec($find, $which);
746 $nearest = $which->dice_axis(1,$found);
747 Returns for each row-vector in $find the index along dimension N of the
749 least row vector of $which greater or equal to it. $which should be
750 sorted in increasing order. If the value of $find is larger than any
751 member of $which, the index to the last element of $which is returned.
752 See also: "vsearch". Contributed by Bryan Jurish <moocow@cpan.org>.
754 vsearchvec does not process bad values. It will set the bad-value flag
756 of all output ndarrays if the flag is set for any of the input
757 ndarrays.
758 unionvec
760 Signature: (a(M,NA); b(M,NB); [o]c(M,NC); indx [o]nc())
761 Union of two vector-valued PDLs.
763 Input PDLs $a() and $b() MUST be sorted in lexicographic order. On
765 return, $nc() holds the actual number of vector-values in the union.
766 In scalar context, slices $c() to the actual number of elements in the
768 union and returns the sliced PDL.
769 Contributed by Bryan Jurish <moocow@cpan.org>.
771 unionvec does not process bad values. It will set the bad-value flag
773 of all output ndarrays if the flag is set for any of the input
774 ndarrays.
775 intersectvec
777 Signature: (a(M,NA); b(M,NB); [o]c(M,NC); indx [o]nc())
778 Intersection of two vector-valued PDLs. Input PDLs $a() and $b() MUST
780 be sorted in lexicographic order. On return, $nc() holds the actual
781 number of vector-values in the intersection.
782 In scalar context, slices $c() to the actual number of elements in the
784 intersection and returns the sliced PDL.
785 Contributed by Bryan Jurish <moocow@cpan.org>.
787 intersectvec does not process bad values. It will set the bad-value
789 flag of all output ndarrays if the flag is set for any of the input
790 ndarrays.
791 setdiffvec
793 Signature: (a(M,NA); b(M,NB); [o]c(M,NC); indx [o]nc())
794 Set-difference ($a() \ $b()) of two vector-valued PDLs.
796 Input PDLs $a() and $b() MUST be sorted in lexicographic order. On
798 return, $nc() holds the actual number of vector-values in the computed
799 vector set.
800 In scalar context, slices $c() to the actual number of elements in the
802 output vector set and returns the sliced PDL.
803 Contributed by Bryan Jurish <moocow@cpan.org>.
805 setdiffvec does not process bad values. It will set the bad-value flag
807 of all output ndarrays if the flag is set for any of the input
808 ndarrays.
809 union_sorted
811 Signature: (a(NA); b(NB); [o]c(NC); indx [o]nc())
812 Union of two flat sorted unique-valued PDLs. Input PDLs $a() and $b()
814 MUST be sorted in lexicographic order and contain no duplicates. On
815 return, $nc() holds the actual number of values in the union.
816 In scalar context, reshapes $c() to the actual number of elements in
818 the union and returns it.
819 Contributed by Bryan Jurish <moocow@cpan.org>.
821 union_sorted does not process bad values. It will set the bad-value
823 flag of all output ndarrays if the flag is set for any of the input
824 ndarrays.
825 intersect_sorted
827 Signature: (a(NA); b(NB); [o]c(NC); indx [o]nc())
828 Intersection of two flat sorted unique-valued PDLs. Input PDLs $a()
830 and $b() MUST be sorted in lexicographic order and contain no
831 duplicates. On return, $nc() holds the actual number of values in the
832 intersection.
833 In scalar context, reshapes $c() to the actual number of elements in
835 the intersection and returns it.
836 Contributed by Bryan Jurish <moocow@cpan.org>.
838 intersect_sorted does not process bad values. It will set the bad-
840 value flag of all output ndarrays if the flag is set for any of the
841 input ndarrays.
842 setdiff_sorted
844 Signature: (a(NA); b(NB); [o]c(NC); indx [o]nc())
845 Set-difference ($a() \ $b()) of two flat sorted unique-valued PDLs.
847 Input PDLs $a() and $b() MUST be sorted in lexicographic order and
849 contain no duplicate values. On return, $nc() holds the actual number
850 of values in the computed vector set.
851 In scalar context, reshapes $c() to the actual number of elements in
853 the difference set and returns it.
854 Contributed by Bryan Jurish <moocow@cpan.org>.
856 setdiff_sorted does not process bad values. It will set the bad-value
858 flag of all output ndarrays if the flag is set for any of the input
859 ndarrays.
860 srand
862 Signature: (a())
863 Seed random-number generator with a 64-bit int. Will generate seed data
865 for a number of threads equal to the return-value of "online_cpus" in
866 PDL::Core.
867 srand(); # uses current time
869 srand(5); # fixed number e.g. for testing
870 srand does not process bad values. It will set the bad-value flag of
872 all output ndarrays if the flag is set for any of the input ndarrays.
873 random
875 Signature: (a())
876 Constructor which returns ndarray of random numbers
878 $x = random([type], $nx, $ny, $nz,...);
880 $x = random $y;
881 etc (see zeroes).
883 This is the uniform distribution between 0 and 1 (assumedly excluding 1
885 itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
886 specify dimensions, types or give a template.
887 You can use the PDL function "srand" to seed the random generator. If
889 it has not been called yet, it will be with the current time.
890 random does not process bad values. It will set the bad-value flag of
892 all output ndarrays if the flag is set for any of the input ndarrays.
893 randsym
895 Signature: (a())
896 Constructor which returns ndarray of random numbers
898 $x = randsym([type], $nx, $ny, $nz,...);
900 $x = randsym $y;
901 etc (see zeroes).
903 This is the uniform distribution between 0 and 1 (excluding both 0 and
905 1, cf "random"). The arguments are the same as "zeroes" (q.v.) - i.e.
906 one can specify dimensions, types or give a template.
907 You can use the PDL function "srand" to seed the random generator. If
909 it has not been called yet, it will be with the current time.
910 randsym does not process bad values. It will set the bad-value flag of
912 all output ndarrays if the flag is set for any of the input ndarrays.
913 grandom
915 Constructor which returns ndarray of Gaussian random numbers
916 $x = grandom([type], $nx, $ny, $nz,...);
918 $x = grandom $y;
919 etc (see zeroes).
921 This is generated using the math library routine "ndtri".
923 Mean = 0, Stddev = 1
925 You can use the PDL function "srand" to seed the random generator. If
927 it has not been called yet, it will be with the current time.
928 vsearch
930 Signature: ( vals(); xs(n); [o] indx(); [\%options] )
931 Efficiently search for values in a sorted ndarray, returning indices.
933 $idx = vsearch( $vals, $x, [\%options] );
935 vsearch( $vals, $x, $idx, [\%options ] );
936 vsearch performs a binary search in the ordered ndarray $x, for the
938 values from $vals ndarray, returning indices into $x. What is a
939 "match", and the meaning of the returned indices, are determined by the
940 options.
941 The "mode" option indicates which method of searching to use, and may
943 be one of:
944 "sample"
946 invoke vsearch_sample, returning indices appropriate for sampling
947 within a distribution.
948 "insert_leftmost"
950 invoke vsearch_insert_leftmost, returning the left-most possible
951 insertion point which still leaves the ndarray sorted.
952 "insert_rightmost"
954 invoke vsearch_insert_rightmost, returning the right-most possible
955 insertion point which still leaves the ndarray sorted.
956 "match"
958 invoke vsearch_match, returning the index of a matching element,
959 else -(insertion point + 1)
960 "bin_inclusive"
962 invoke vsearch_bin_inclusive, returning an index appropriate for
963 binning on a grid where the left bin edges are inclusive of the
964 bin. See below for further explanation of the bin.
965 "bin_exclusive"
967 invoke vsearch_bin_exclusive, returning an index appropriate for
968 binning on a grid where the left bin edges are exclusive of the
969 bin. See below for further explanation of the bin.
970 The default value of "mode" is "sample".
972 use PDL;
974 my @modes = qw( sample insert_leftmost insert_rightmost match
976 bin_inclusive bin_exclusive );
977 # Generate a sequence of 3 zeros, 3 ones, ..., 3 fours.
979 my $x = zeroes(3,5)->yvals->flat;
980 for my $mode ( @modes ) {
982 # if the value is in $x
983 my $contained = 2;
984 my $idx_contained = vsearch( $contained, $x, { mode => $mode } );
985 my $x_contained = $x->copy;
986 $x_contained->slice( $idx_contained ) .= 9;
987 # if the value is not in $x
989 my $not_contained = 1.5;
990 my $idx_not_contained = vsearch( $not_contained, $x, { mode => $mode } );
991 my $x_not_contained = $x->copy;
992 $x_not_contained->slice( $idx_not_contained ) .= 9;
993 print sprintf("%-23s%30s\n", '$x', $x);
995 print sprintf("%-23s%30s\n", "$mode ($contained)", $x_contained);
996 print sprintf("%-23s%30s\n\n", "$mode ($not_contained)", $x_not_contained);
997 }
998 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
1000 # sample (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
1001 # sample (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
1002 #
1003 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
1004 # insert_leftmost (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
1005 # insert_leftmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
1006 #
1007 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
1008 # insert_rightmost (2) [0 0 0 1 1 1 2 2 2 9 3 3 4 4 4]
1009 # insert_rightmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
1010 #
1011 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
1012 # match (2) [0 0 0 1 1 1 2 9 2 3 3 3 4 4 4]
1013 # match (1.5) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
1014 #
1015 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
1016 # bin_inclusive (2) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
1017 # bin_inclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
1018 #
1019 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
1020 # bin_exclusive (2) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
1021 # bin_exclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
1022 Also see vsearch_sample, vsearch_insert_leftmost,
1024 vsearch_insert_rightmost, vsearch_match, vsearch_bin_inclusive, and
1025 vsearch_bin_exclusive
1026 vsearch_sample
1028 Signature: (vals(); x(n); indx [o]idx())
1029 Search for values in a sorted array, return index appropriate for
1031 sampling from a distribution
1032 $idx = vsearch_sample($vals, $x);
1034 $x must be sorted, but may be in decreasing or increasing order.
1036 vsearch_sample returns an index I for each value V of $vals appropriate
1038 for sampling $vals
1039 I has the following properties:
1041 • if $x is sorted in increasing order
1043 V <= x[0] : I = 0
1045 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
1046 x[-1] < V : I = $x->nelem -1
1047 • if $x is sorted in decreasing order
1049 V > x[0] : I = 0
1051 x[0] >= V > x[-1] : I s.t. x[I] >= V > x[I+1]
1052 x[-1] >= V : I = $x->nelem - 1
1053 If all elements of $x are equal, I = $x->nelem - 1.
1055 If $x contains duplicated elements, I is the index of the leftmost (by
1057 position in array) duplicate if V matches.
1058 This function is useful e.g. when you have a list of probabilities for
1060 events and want to generate indices to events:
1061 $x = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
1063 $y = random 20;
1064 $c = vsearch_sample($y, $x); # Now, $c will have the appropriate distr.
1065 It is possible to use the cumusumover function to obtain cumulative
1067 probabilities from absolute probabilities.
1068 needs major (?) work to handles bad values
1070 vsearch_insert_leftmost
1072 Signature: (vals(); x(n); indx [o]idx())
1073 Determine the insertion point for values in a sorted array, inserting
1075 before duplicates.
1076 $idx = vsearch_insert_leftmost($vals, $x);
1078 $x must be sorted, but may be in decreasing or increasing order.
1080 vsearch_insert_leftmost returns an index I for each value V of $vals
1082 equal to the leftmost position (by index in array) within $x that V may
1083 be inserted and still maintain the order in $x.
1084 Insertion at index I involves shifting elements I and higher of $x to
1086 the right by one and setting the now empty element at index I to V.
1087 I has the following properties:
1089 • if $x is sorted in increasing order
1091 V <= x[0] : I = 0
1093 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
1094 x[-1] < V : I = $x->nelem
1095 • if $x is sorted in decreasing order
1097 V > x[0] : I = -1
1099 x[0] >= V >= x[-1] : I s.t. x[I] >= V > x[I+1]
1100 x[-1] >= V : I = $x->nelem -1
1101 If all elements of $x are equal,
1103 i = 0
1105 If $x contains duplicated elements, I is the index of the leftmost (by
1107 index in array) duplicate if V matches.
1108 needs major (?) work to handles bad values
1110 vsearch_insert_rightmost
1112 Signature: (vals(); x(n); indx [o]idx())
1113 Determine the insertion point for values in a sorted array, inserting
1115 after duplicates.
1116 $idx = vsearch_insert_rightmost($vals, $x);
1118 $x must be sorted, but may be in decreasing or increasing order.
1120 vsearch_insert_rightmost returns an index I for each value V of $vals
1122 equal to the rightmost position (by index in array) within $x that V
1123 may be inserted and still maintain the order in $x.
1124 Insertion at index I involves shifting elements I and higher of $x to
1126 the right by one and setting the now empty element at index I to V.
1127 I has the following properties:
1129 • if $x is sorted in increasing order
1131 V < x[0] : I = 0
1133 x[0] <= V < x[-1] : I s.t. x[I-1] <= V < x[I]
1134 x[-1] <= V : I = $x->nelem
1135 • if $x is sorted in decreasing order
1137 V >= x[0] : I = -1
1139 x[0] > V >= x[-1] : I s.t. x[I] >= V > x[I+1]
1140 x[-1] > V : I = $x->nelem -1
1141 If all elements of $x are equal,
1143 i = $x->nelem - 1
1145 If $x contains duplicated elements, I is the index of the leftmost (by
1147 index in array) duplicate if V matches.
1148 needs major (?) work to handles bad values
1150 vsearch_match
1152 Signature: (vals(); x(n); indx [o]idx())
1153 Match values against a sorted array.
1155 $idx = vsearch_match($vals, $x);
1157 $x must be sorted, but may be in decreasing or increasing order.
1159 vsearch_match returns an index I for each value V of $vals. If V
1161 matches an element in $x, I is the index of that element, otherwise it
1162 is -( insertion_point + 1 ), where insertion_point is an index in $x
1163 where V may be inserted while maintaining the order in $x. If $x has
1164 duplicated values, I may refer to any of them.
1165 needs major (?) work to handles bad values
1167 vsearch_bin_inclusive
1169 Signature: (vals(); x(n); indx [o]idx())
1170 Determine the index for values in a sorted array of bins, lower bound
1172 inclusive.
1173 $idx = vsearch_bin_inclusive($vals, $x);
1175 $x must be sorted, but may be in decreasing or increasing order.
1177 $x represents the edges of contiguous bins, with the first and last
1179 elements representing the outer edges of the outer bins, and the inner
1180 elements the shared bin edges.
1181 The lower bound of a bin is inclusive to the bin, its outer bound is
1183 exclusive to it. vsearch_bin_inclusive returns an index I for each
1184 value V of $vals
1185 I has the following properties:
1187 • if $x is sorted in increasing order
1189 V < x[0] : I = -1
1191 x[0] <= V < x[-1] : I s.t. x[I] <= V < x[I+1]
1192 x[-1] <= V : I = $x->nelem - 1
1193 • if $x is sorted in decreasing order
1195 V >= x[0] : I = 0
1197 x[0] > V >= x[-1] : I s.t. x[I+1] > V >= x[I]
1198 x[-1] > V : I = $x->nelem
1199 If all elements of $x are equal,
1201 i = $x->nelem - 1
1203 If $x contains duplicated elements, I is the index of the righmost (by
1205 index in array) duplicate if V matches.
1206 needs major (?) work to handles bad values
1208 vsearch_bin_exclusive
1210 Signature: (vals(); x(n); indx [o]idx())
1211 Determine the index for values in a sorted array of bins, lower bound
1213 exclusive.
1214 $idx = vsearch_bin_exclusive($vals, $x);
1216 $x must be sorted, but may be in decreasing or increasing order.
1218 $x represents the edges of contiguous bins, with the first and last
1220 elements representing the outer edges of the outer bins, and the inner
1221 elements the shared bin edges.
1222 The lower bound of a bin is exclusive to the bin, its upper bound is
1224 inclusive to it. vsearch_bin_exclusive returns an index I for each
1225 value V of $vals.
1226 I has the following properties:
1228 • if $x is sorted in increasing order
1230 V <= x[0] : I = -1
1232 x[0] < V <= x[-1] : I s.t. x[I] < V <= x[I+1]
1233 x[-1] < V : I = $x->nelem - 1
1234 • if $x is sorted in decreasing order
1236 V > x[0] : I = 0
1238 x[0] >= V > x[-1] : I s.t. x[I-1] >= V > x[I]
1239 x[-1] >= V : I = $x->nelem
1240 If all elements of $x are equal,
1242 i = $x->nelem - 1
1244 If $x contains duplicated elements, I is the index of the righmost (by
1246 index in array) duplicate if V matches.
1247 needs major (?) work to handles bad values
1249 interpolate
1251 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
1252 routine for 1D linear interpolation
1254 ( $yi, $err ) = interpolate($xi, $x, $y)
1256 Given a set of points "($x,$y)", use linear interpolation to find the
1258 values $yi at a set of points $xi.
1259 "interpolate" uses a binary search to find the suspects, er...,
1261 interpolation indices and therefore abscissas (ie $x) have to be
1262 strictly ordered (increasing or decreasing). For interpolation at lots
1263 of closely spaced abscissas an approach that uses the last index found
1264 as a start for the next search can be faster (compare Numerical Recipes
1265 "hunt" routine). Feel free to implement that on top of the binary
1266 search if you like. For out of bounds values it just does a linear
1267 extrapolation and sets the corresponding element of $err to 1, which is
1268 otherwise 0.
1269 See also "interpol", which uses the same routine, differing only in the
1271 handling of extrapolation - an error message is printed rather than
1272 returning an error ndarray.
1273 needs major (?) work to handles bad values
1275 interpol
1277 Signature: (xi(); x(n); y(n); [o] yi())
1278 routine for 1D linear interpolation
1280 $yi = interpol($xi, $x, $y)
1282 "interpol" uses the same search method as "interpolate", hence $x must
1284 be strictly ordered (either increasing or decreasing). The difference
1285 occurs in the handling of out-of-bounds values; here an error message
1286 is printed.
1287 interpND
1289 Interpolate values from an N-D ndarray, with switchable method
1290 $source = 10*xvals(10,10) + yvals(10,10);
1292 $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
1293 print $source->interpND( $index );
1294 InterpND acts like indexND, collapsing $index by lookup into $source;
1296 but it does interpolation rather than direct sampling. The
1297 interpolation method and boundary condition are switchable via an
1298 options hash.
1299 By default, linear or sample interpolation is used, with constant value
1301 outside the boundaries of the source pdl. No dataflow occurs, because
1302 in general the output is computed rather than indexed.
1303 All the interpolation methods treat the pixels as value-centered, so
1305 the "sample" method will return "$a->(0)" for coordinate values on the
1306 set [-0.5,0.5), and all methods will return "$a->(1)" for a coordinate
1307 value of exactly 1.
1308 Recognized options:
1310 method
1312 Values can be:
1313 • 0, s, sample, Sample (default for integer source types)
1315 The nearest value is taken. Pixels are regarded as centered on
1317 their respective integer coordinates (no offset from the linear
1318 case).
1319 • 1, l, linear, Linear (default for floating point source types)
1321 The values are N-linearly interpolated from an N-dimensional cube
1323 of size 2.
1324 • 3, c, cube, cubic, Cubic
1326 The values are interpolated using a local cubic fit to the data.
1328 The fit is constrained to match the original data and its
1329 derivative at the data points. The second derivative of the fit
1330 is not continuous at the data points. Multidimensional datasets
1331 are interpolated by the successive-collapse method.
1332 (Note that the constraint on the first derivative causes a small
1334 amount of ringing around sudden features such as step functions).
1335 • f, fft, fourier, Fourier
1337 The source is Fourier transformed, and the interpolated values
1339 are explicitly calculated from the coefficients. The boundary
1340 condition option is ignored -- periodic boundaries are imposed.
1341 If you pass in the option "fft", and it is a list (ARRAY) ref,
1343 then it is a stash for the magnitude and phase of the source FFT.
1344 If the list has two elements then they are taken as already
1345 computed; otherwise they are calculated and put in the stash.
1346 b, bound, boundary, Boundary
1348 This option is passed unmodified into indexND, which is used as the
1349 indexing engine for the interpolation. Some current allowed values
1350 are 'extend', 'periodic', 'truncate', and 'mirror' (default is
1351 'truncate').
1352 bad
1354 contains the fill value used for 'truncate' boundary. (default 0)
1355 fft
1357 An array ref whose associated list is used to stash the FFT of the
1358 source data, for the FFT method.
1359 one2nd
1361 Converts a one dimensional index ndarray to a set of ND coordinates
1362 @coords=one2nd($x, $indices)
1364 returns an array of ndarrays containing the ND indexes corresponding to
1366 the one dimensional list indices. The indices are assumed to correspond
1367 to array $x clumped using "clump(-1)". This routine is used in the old
1368 vector form of "whichND", but is useful on its own occasionally.
1369 Returned ndarrays have the indx datatype. $indices can have values
1371 larger than "$x->nelem" but negative values in $indices will not give
1372 the answer you expect.
1373 pdl> $x=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$x->clump(-1)
1375 pdl> $maxind=maximum_ind($c); p $maxind;
1376 6
1377 pdl> print one2nd($x, maximum_ind($c))
1378 0 1 1
1379 pdl> p $x->at(0,1,1)
1380 3
1381 which
1383 Signature: (mask(n); indx [o] inds(n); indx [o]lastout())
1384 Returns indices of non-zero values from a 1-D PDL
1386 $i = which($mask);
1388 returns a pdl with indices for all those elements that are nonzero in
1390 the mask. Note that the returned indices will be 1D. If you feed in a
1391 multidimensional mask, it will be flattened before the indices are
1392 calculated. See also "whichND" for multidimensional masks.
1393 If you want to index into the original mask or a similar ndarray with
1395 output from "which", remember to flatten it before calling index:
1396 $data = random 5, 5;
1398 $idx = which $data > 0.5; # $idx is now 1D
1399 $bigsum = $data->flat->index($idx)->sum; # flatten before indexing
1400 Compare also "where" for similar functionality.
1402 SEE ALSO:
1404 "which_both" returns separately the indices of both nonzero and zero
1406 values in the mask.
1407 "where_both" returns separately slices of both nonzero and zero values
1409 in the mask.
1410 "where" returns associated values from a data PDL, rather than indices
1412 into the mask PDL.
1413 "whichND" returns N-D indices into a multidimensional PDL.
1415 pdl> $x = sequence(10); p $x
1417 [0 1 2 3 4 5 6 7 8 9]
1418 pdl> $indx = which($x>6); p $indx
1419 [7 8 9]
1420 which processes bad values. It will set the bad-value flag of all
1422 output ndarrays if the flag is set for any of the input ndarrays.
1423 which_both
1425 Signature: (mask(n); indx [o] inds(n); indx [o]notinds(n); indx [o]lastout(); indx [o]lastoutn())
1426 Returns indices of nonzero and zero values in a mask PDL
1428 ($i, $c_i) = which_both($mask);
1430 This works just as "which", but the complement of $i will be in $c_i.
1432 pdl> p $x = sequence(10)
1434 [0 1 2 3 4 5 6 7 8 9]
1435 pdl> ($big, $small) = which_both($x >= 5); p "$big\n$small"
1436 [5 6 7 8 9]
1437 [0 1 2 3 4]
1438 which_both processes bad values. It will set the bad-value flag of all
1440 output ndarrays if the flag is set for any of the input ndarrays.
1441 where
1443 Use a mask to select values from one or more data PDLs
1444 "where" accepts one or more data ndarrays and a mask ndarray. It
1446 returns a list of output ndarrays, corresponding to the input data
1447 ndarrays. Each output ndarray is a 1-dimensional list of values in its
1448 corresponding data ndarray. The values are drawn from locations where
1449 the mask is nonzero.
1450 The output PDLs are still connected to the original data PDLs, for the
1452 purpose of dataflow.
1453 "where" combines the functionality of "which" and index into a single
1455 operation.
1456 BUGS:
1458 While "where" works OK for most N-dimensional cases, it does not
1460 broadcast properly over (for example) the (N+1)th dimension in data
1461 that is compared to an N-dimensional mask. Use "whereND" for that.
1462 $i = $x->where($x+5 > 0); # $i contains those elements of $x
1464 # where mask ($x+5 > 0) is 1
1465 $i .= -5; # Set those elements (of $x) to -5. Together, these
1466 # commands clamp $x to a maximum of -5.
1467 It is also possible to use the same mask for several ndarrays with the
1469 same call:
1470 ($i,$j,$k) = where($x,$y,$z, $x+5>0);
1472 Note: $i is always 1-D, even if $x is >1-D.
1474 WARNING: The first argument (the values) and the second argument (the
1476 mask) currently have to have the exact same dimensions (or horrible
1477 things happen). You *cannot* broadcast over a smaller mask, for
1478 example.
1479 where_both
1481 Returns slices (non-zero in mask, zero) of an ndarray according to a
1482 mask
1483 ($match_vals, $non_match_vals) = where_both($pdl, $mask);
1485 This works like "which_both", but (flattened) data-flowing slices
1487 rather than index-sets are returned.
1488 pdl> p $x = sequence(10) + 2
1490 [2 3 4 5 6 7 8 9 10 11]
1491 pdl> ($big, $small) = where_both($x, $x > 5); p "$big\n$small"
1492 [6 7 8 9 10 11]
1493 [2 3 4 5]
1494 pdl> p $big += 2, $small -= 1
1495 [8 9 10 11 12 13] [1 2 3 4]
1496 pdl> p $x
1497 [1 2 3 4 8 9 10 11 12 13]
1498 whereND
1500 "where" with support for ND masks and broadcasting
1501 "whereND" accepts one or more data ndarrays and a mask ndarray. It
1503 returns a list of output ndarrays, corresponding to the input data
1504 ndarrays. The values are drawn from locations where the mask is
1505 nonzero.
1506 "whereND" differs from "where" in that the mask dimensionality is
1508 preserved which allows for proper broadcasting of the selection
1509 operation over higher dimensions.
1510 As with "where" the output PDLs are still connected to the original
1512 data PDLs, for the purpose of dataflow.
1513 $sdata = whereND $data, $mask
1515 ($s1, $s2, ..., $sn) = whereND $d1, $d2, ..., $dn, $mask
1516 where
1518 $data is M dimensional
1520 $mask is N < M dimensional
1521 dims($data) 1..N == dims($mask) 1..N
1522 with broadcasting over N+1 to M dimensions
1523 $data = sequence(4,3,2); # example data array
1525 $mask4 = (random(4)>0.5); # example 1-D mask array, has $n4 true values
1526 $mask43 = (random(4,3)>0.5); # example 2-D mask array, has $n43 true values
1527 $sdat4 = whereND $data, $mask4; # $sdat4 is a [$n4,3,2] pdl
1528 $sdat43 = whereND $data, $mask43; # $sdat43 is a [$n43,2] pdl
1529 Just as with "where", you can use the returned value in an assignment.
1531 That means that both of these examples are valid:
1532 # Used to create a new slice stored in $sdat4:
1534 $sdat4 = $data->whereND($mask4);
1535 $sdat4 .= 0;
1536 # Used in lvalue context:
1537 $data->whereND($mask4) .= 0;
1538 SEE ALSO:
1540 "whichND" returns N-D indices into a multidimensional PDL, from a mask.
1542 whichND
1544 Return the coordinates of non-zero values in a mask.
1545 WhichND returns the N-dimensional coordinates of each nonzero value in
1547 a mask PDL with any number of dimensions. The returned values arrive
1548 as an array-of-vectors suitable for use in indexND or range.
1549 $coords = whichND($mask);
1551 returns a PDL containing the coordinates of the elements that are non-
1553 zero in $mask, suitable for use in "indexND" in PDL::Slices. The 0th
1554 dimension contains the full coordinate listing of each point; the 1st
1555 dimension lists all the points. For example, if $mask has rank 4 and
1556 100 matching elements, then $coords has dimension 4x100.
1557 If no such elements exist, then whichND returns a structured empty PDL:
1559 an Nx0 PDL that contains no values (but matches, broadcasting-wise,
1560 with the vectors that would be produced if such elements existed).
1561 DEPRECATED BEHAVIOR IN LIST CONTEXT:
1563 whichND once delivered different values in list context than in scalar
1565 context, for historical reasons. In list context, it returned the
1566 coordinates transposed, as a collection of 1-PDLs (one per dimension)
1567 in a list. This usage is deprecated in PDL 2.4.10, and will cause a
1568 warning to be issued every time it is encountered. To avoid the
1569 warning, you can set the global variable "$PDL::whichND" to 's' to get
1570 scalar behavior in all contexts, or to 'l' to get list behavior in list
1571 context.
1572 In later versions of PDL, the deprecated behavior will disappear.
1574 Deprecated list context whichND expressions can be replaced with:
1575 @list = $x->whichND->mv(0,-1)->dog;
1577 SEE ALSO:
1579 "which" finds coordinates of nonzero values in a 1-D mask.
1581 "where" extracts values from a data PDL that are associated with
1583 nonzero values in a mask PDL.
1584 "indexND" in PDL::Slices can be fed the coordinates to return the
1586 values.
1587 pdl> $s=sequence(10,10,3,4)
1589 pdl> ($x, $y, $z, $w)=whichND($s == 203); p $x, $y, $z, $w
1590 [3] [0] [2] [0]
1591 pdl> print $s->at(list(cat($x,$y,$z,$w)))
1592 203
1593 setops
1595 Implements simple set operations like union and intersection
1596 Usage: $set = setops($x, <OPERATOR>, $y);
1598 The operator can be "OR", "XOR" or "AND". This is then applied to $x
1600 viewed as a set and $y viewed as a set. Set theory says that a set may
1601 not have two or more identical elements, but setops takes care of this
1602 for you, so "$x=pdl(1,1,2)" is OK. The functioning is as follows:
1603 "OR"
1605 The resulting vector will contain the elements that are either in
1606 $x or in $y or both. This is the union in set operation terms
1607 "XOR"
1609 The resulting vector will contain the elements that are either in
1610 $x or $y, but not in both. This is
1611 Union($x, $y) - Intersection($x, $y)
1613 in set operation terms.
1615 "AND"
1617 The resulting vector will contain the intersection of $x and $y, so
1618 the elements that are in both $x and $y. Note that for convenience
1619 this operation is also aliased to "intersect".
1620 It should be emphasized that these routines are used when one or both
1622 of the sets $x, $y are hard to calculate or that you get from a
1623 separate subroutine.
1624 Finally IDL users might be familiar with Craig Markwardt's
1626 "cmset_op.pro" routine which has inspired this routine although it was
1627 written independently However the present routine has a few less
1628 options (but see the examples)
1629 You will very often use these functions on an index vector, so that is
1631 what we will show here. We will in fact something slightly silly. First
1632 we will find all squares that are also cubes below 10000.
1633 Create a sequence vector:
1635 pdl> $x = sequence(10000)
1637 Find all odd and even elements:
1639 pdl> ($even, $odd) = which_both( ($x % 2) == 0)
1641 Find all squares
1643 pdl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
1645 Find all cubes (being careful with roundoff error!)
1647 pdl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
1649 Then find all squares that are cubes:
1651 pdl> $both = setops($squares, 'AND', $cubes)
1653 And print these (assumes that "PDL::NiceSlice" is loaded!)
1655 pdl> p $x($both)
1657 [0 1 64 729 4096]
1658 Then find all numbers that are either cubes or squares, but not both:
1660 pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)
1662 pdl> p $cube_xor_square->nelem()
1664 112
1665 So there are a total of 112 of these!
1667 Finally find all odd squares:
1669 pdl> $odd_squares = setops($squares, 'AND', $odd)
1671 Another common occurrence is to want to get all objects that are in $x
1673 and in the complement of $y. But it is almost always best to create the
1674 complement explicitly since the universe that both are taken from is
1675 not known. Thus use "which_both" if possible to keep track of
1676 complements.
1677 If this is impossible the best approach is to make a temporary:
1679 This creates an index vector the size of the universe of the sets and
1681 set all elements in $y to 0
1682 pdl> $tmp = ones($n_universe); $tmp($y) .= 0;
1684 This then finds the complement of $y
1686 pdl> $C_b = which($tmp == 1);
1688 and this does the final selection:
1690 pdl> $set = setops($x, 'AND', $C_b)
1692 intersect
1694 Calculate the intersection of two ndarrays
1695 Usage: $set = intersect($x, $y);
1697 This routine is merely a simple interface to "setops". See that for
1699 more information
1700 Find all numbers less that 100 that are of the form 2*y and 3*x
1702 pdl> $x=sequence(100)
1704 pdl> $factor2 = which( ($x % 2) == 0)
1705 pdl> $factor3 = which( ($x % 3) == 0)
1706 pdl> $ii=intersect($factor2, $factor3)
1707 pdl> p $x($ii)
1708 [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
1709
1711 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
1712 Contributions by Christian Soeller (c.soeller@auckland.ac.nz), Karl
1713 Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
1714 (deforest@boulder.swri.edu) and Jarle Brinchmann (jarle@astro.up.pt)
1715 All rights reserved. There is no warranty. You are allowed to
1716 redistribute this software / documentation under certain conditions.
1717 For details, see the file COPYING in the PDL distribution. If this file
1718 is separated from the PDL distribution, the copyright notice should be
1719 included in the file.
1720 Updated for CPAN viewing compatibility by David Mertens.
1722perl v5.36.0 2022-07-22 Primitive(3)