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: (a(n); indx ind(n); [o] sum(m))
229 Threaded Index Add: Add "a" to the "ind" element of "sum", i.e:
231 sum(ind) += a
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 Threaded 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 if any of the indices are bad.
255 conv1d
257 Signature: (a(m); kern(p); [o]b(m); int reflect)
258 1D convolution along first dimension
260 The m-th element of the discrete convolution of an input ndarray $a of
262 size $M, and a kernel ndarray $kern of size $P, is calculated as
263 n = ($P-1)/2
265 ====
266 \
267 ($a conv1d $kern)[m] = > $a_ext[m - n] * $kern[n]
268 /
269 ====
270 n = -($P-1)/2
271 where $a_ext is either the periodic (or reflected) extension of $a so
273 it is equal to $a on " 0..$M-1 " and equal to the corresponding
274 periodic/reflected image of $a outside that range.
275 $con = conv1d sequence(10), pdl(-1,0,1);
277 $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
279 By default, periodic boundary conditions are assumed (i.e. wrap
281 around). Alternatively, you can request reflective boundary conditions
282 using the "Boundary" option:
283 {Boundary => 'reflect'} # case in 'reflect' doesn't matter
285 The convolution is performed along the first dimension. To apply it
287 across another dimension use the slicing routines, e.g.
288 $y = $x->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
290 This function is useful for broadcasted filtering of 1D signals.
292 Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv
294 WARNING: "conv1d" processes bad values in its inputs as the numeric
296 value of "$pdl->badvalue" so it is not recommended for processing pdls
297 with bad values in them unless special care is taken.
298 conv1d ignores the bad-value flag of the input ndarrays. It will set
300 the bad-value flag of all output ndarrays if the flag is set for any of
301 the input ndarrays.
302 in
304 Signature: (a(); b(n); [o] c())
305 test if a is in the set of values b
307 $goodmsk = $labels->in($goodlabels);
309 print pdl(3,1,4,6,2)->in(pdl(2,3,3));
310 [1 0 0 0 1]
311 "in" is akin to the is an element of of set theory. In principle, PDL
313 broadcasting could be used to achieve its functionality by using a
314 construct like
315 $msk = ($labels->dummy(0) == $goodlabels)->orover;
317 However, "in" doesn't create a (potentially large) intermediate and is
319 generally faster.
320 in does not process bad values. It will set the bad-value flag of all
322 output ndarrays if the flag is set for any of the input ndarrays.
323 uniq
325 return all unique elements of an ndarray
326 The unique elements are returned in ascending order.
328 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniq
330 [-1 0 2 4 6] # 0 is returned 2nd (sorted order)
331 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniq
333 [-1 2 4 6 nan] # NaN value is returned at end
334 Note: The returned pdl is 1D; any structure of the input ndarray is
336 lost. "NaN" values are never compare equal to any other values, even
337 themselves. As a result, they are always unique. "uniq" returns the
338 NaN values at the end of the result ndarray. This follows the Matlab
339 usage.
340 See "uniqind" if you need the indices of the unique elements rather
342 than the values.
343 Bad values are not considered unique by uniq and are ignored.
345 $x=sequence(10);
347 $x=$x->setbadif($x%3);
348 print $x->uniq;
349 [0 3 6 9]
350 uniqind
352 Return the indices of all unique elements of an ndarray The order is in
353 the order of the values to be consistent with uniq. "NaN" values never
354 compare equal with any other value and so are always unique. This
355 follows the Matlab usage.
356 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniqind
358 [5 4 1 3 6] # the 0 at index 4 is returned 2nd, but...
359 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniqind
361 [5 1 3 6 4] # ...the NaN at index 4 is returned at end
362 Note: The returned pdl is 1D; any structure of the input ndarray is
364 lost.
365 See "uniq" if you want the unique values instead of the indices.
367 Bad values are not considered unique by uniqind and are ignored.
369 uniqvec
371 Return all unique vectors out of a collection
372 NOTE: If any vectors in the input ndarray have NaN values
374 they are returned at the end of the non-NaN ones. This is
375 because, by definition, NaN values never compare equal with
376 any other value.
377 NOTE: The current implementation does not sort the vectors
379 containing NaN values.
380 The unique vectors are returned in lexicographically sorted ascending
382 order. The 0th dimension of the input PDL is treated as a dimensional
383 index within each vector, and the 1st and any higher dimensions are
384 taken to run across vectors. The return value is always 2D; any
385 structure of the input PDL (beyond using the 0th dimension for vector
386 index) is lost.
387 See also "uniq" for a unique list of scalars; and qsortvec for sorting
389 a list of vectors lexicographcally.
390 If a vector contains all bad values, it is ignored as in "uniq". If
392 some of the values are good, it is treated as a normal vector. For
393 example, [1 2 BAD] and [BAD 2 3] could be returned, but [BAD BAD BAD]
394 could not. Vectors containing BAD values will be returned after any
395 non-NaN and non-BAD containing vectors, followed by the NaN vectors.
396 hclip
398 Signature: (a(); b(); [o] c())
399 clip (threshold) $a by $b ($b is upper bound)
401 hclip processes bad values. It will set the bad-value flag of all
403 output ndarrays if the flag is set for any of the input ndarrays.
404 lclip
406 Signature: (a(); b(); [o] c())
407 clip (threshold) $a by $b ($b is lower bound)
409 lclip processes bad values. It will set the bad-value flag of all
411 output ndarrays if the flag is set for any of the input ndarrays.
412 clip
414 Clip (threshold) an ndarray by (optional) upper or lower bounds.
415 $y = $x->clip(0,3);
417 $c = $x->clip(undef, $x);
418 clip handles bad values since it is just a wrapper around "hclip" and
420 "lclip".
421 clip
423 Signature: (a(); l(); h(); [o] c())
424 info not available
426 clip processes bad values. It will set the bad-value flag of all
428 output ndarrays if the flag is set for any of the input ndarrays.
429 wtstat
431 Signature: (a(n); wt(n); avg(); [o]b(); int deg)
432 Weighted statistical moment of given degree
434 This calculates a weighted statistic over the vector "a". The formula
436 is
437 b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
439 Bad values are ignored in any calculation; $b will only have its bad
441 flag set if the output contains any bad data.
442 statsover
444 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())
445 Calculate useful statistics over a dimension of an ndarray
447 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($ndarray, $weights);
449 This utility function calculates various useful quantities of an
451 ndarray. These are:
452 • the mean:
454 MEAN = sum (x)/ N
456 with "N" being the number of elements in x
458 • the population RMS deviation from the mean:
460 PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
462 The population deviation is the best-estimate of the deviation of
464 the population from which a sample is drawn.
465 • the median
467 The median is the 50th percentile data value. Median is found by
469 medover, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
470 • the minimum
472 • the maximum
474 • the average absolute deviation:
476 AADEV = sum( abs(x-mean(x)) )/N
478 • RMS deviation from the mean:
480 RMS = sqrt(sum( (x-mean(x))^2 )/N)
482 (also known as the root-mean-square deviation, or the square root of
484 the variance)
485 This operator is a projection operator so the calculation will take
487 place over the final dimension. Thus if the input is N-dimensional each
488 returned value will be N-1 dimensional, to calculate the statistics for
489 the entire ndarray either use "clump(-1)" directly on the ndarray or
490 call "stats".
491 Bad values are simply ignored in the calculation, effectively reducing
493 the sample size. If all data are bad then the output data are marked
494 bad.
495 stats
497 Calculates useful statistics on an ndarray
498 ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($ndarray,[$weights]);
500 This utility calculates all the most useful quantities in one call. It
502 works the same way as "statsover", except that the quantities are
503 calculated considering the entire input PDL as a single sample, rather
504 than as a collection of rows. See "statsover" for definitions of the
505 returned quantities.
506 Bad values are handled; if all input values are bad, then all of the
508 output values are flagged bad.
509 histogram
511 Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
512 Calculates a histogram for given stepsize and minimum.
514 $h = histogram($data, $step, $min, $numbins);
516 $hist = zeroes $numbins; # Put histogram in existing ndarray.
517 histogram($data, $hist, $step, $min, $numbins);
518 The histogram will contain $numbins bins starting from $min, each $step
520 wide. The value in each bin is the number of values in $data that lie
521 within the bin limits.
522 Data below the lower limit is put in the first bin, and data above the
524 upper limit is put in the last bin.
525 The output is reset in a different broadcastloop so that you can take a
527 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
528 For a higher-level interface, see hist.
530 pdl> p histogram(pdl(1,1,2),1,0,3)
532 [0 2 1]
533 histogram processes bad values. It will set the bad-value flag of all
535 output ndarrays if the flag is set for any of the input ndarrays.
536 whistogram
538 Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
539 Calculates a histogram from weighted data for given stepsize and
541 minimum.
542 $h = whistogram($data, $weights, $step, $min, $numbins);
544 $hist = zeroes $numbins; # Put histogram in existing ndarray.
545 whistogram($data, $weights, $hist, $step, $min, $numbins);
546 The histogram will contain $numbins bins starting from $min, each $step
548 wide. The value in each bin is the sum of the values in $weights that
549 correspond to values in $data that lie within the bin limits.
550 Data below the lower limit is put in the first bin, and data above the
552 upper limit is put in the last bin.
553 The output is reset in a different broadcastloop so that you can take a
555 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
556 pdl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
558 [0 0.2 0.5 0]
559 whistogram processes bad values. It will set the bad-value flag of all
561 output ndarrays if the flag is set for any of the input ndarrays.
562 histogram2d
564 Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
565 double stepb; double minb; int mbsize => mb;)
566 Calculates a 2d histogram.
568 $h = histogram2d($datax, $datay, $stepx, $minx,
570 $nbinx, $stepy, $miny, $nbiny);
571 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing ndarray.
572 histogram2d($datax, $datay, $hist, $stepx, $minx,
573 $nbinx, $stepy, $miny, $nbiny);
574 The histogram will contain $nbinx x $nbiny bins, with the lower limits
576 of the first one at "($minx, $miny)", and with bin size "($stepx,
577 $stepy)". The value in each bin is the number of values in $datax and
578 $datay that lie within the bin limits.
579 Data below the lower limit is put in the first bin, and data above the
581 upper limit is put in the last bin.
582 pdl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
584 [
585 [0 0 0]
586 [0 2 2]
587 [0 1 0]
588 ]
589 histogram2d processes bad values. It will set the bad-value flag of
591 all output ndarrays if the flag is set for any of the input ndarrays.
592 whistogram2d
594 Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
595 double stepb; double minb; int mbsize => mb;)
596 Calculates a 2d histogram from weighted data.
598 $h = whistogram2d($datax, $datay, $weights,
600 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
601 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing ndarray.
602 whistogram2d($datax, $datay, $weights, $hist,
603 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
604 The histogram will contain $nbinx x $nbiny bins, with the lower limits
606 of the first one at "($minx, $miny)", and with bin size "($stepx,
607 $stepy)". The value in each bin is the sum of the values in $weights
608 that correspond to values in $datax and $datay that lie within the bin
609 limits.
610 Data below the lower limit is put in the first bin, and data above the
612 upper limit is put in the last bin.
613 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)
615 [
616 [ 0 0 0]
617 [ 0 0.5 0.9]
618 [ 0 0.1 0]
619 ]
620 whistogram2d processes bad values. It will set the bad-value flag of
622 all output ndarrays if the flag is set for any of the input ndarrays.
623 fibonacci
625 Signature: (i(n); indx [o]x(n))
626 Constructor - a vector with Fibonacci's sequence
628 fibonacci does not process bad values. It will set the bad-value flag
630 of all output ndarrays if the flag is set for any of the input
631 ndarrays.
632 append
634 Signature: (a(n); b(m); [o] c(mn))
635 append two ndarrays by concatenating along their first dimensions
637 $x = ones(2,4,7);
639 $y = sequence 5;
640 $c = $x->append($y); # size of $c is now (7,4,7) (a jumbo-ndarray ;)
641 "append" appends two ndarrays along their first dimensions. The rest of
643 the dimensions must be compatible in the broadcasting sense. The
644 resulting size of the first dimension is the sum of the sizes of the
645 first dimensions of the two argument ndarrays - i.e. "n + m".
646 Similar functions include "glue" (below), which can append more than
648 two ndarrays along an arbitrary dimension, and cat, which can append
649 more than two ndarrays that all have the same sized dimensions.
650 append does not process bad values. It will set the bad-value flag of
652 all output ndarrays if the flag is set for any of the input ndarrays.
653 glue
655 $c = $x->glue(<dim>,$y,...)
656 Glue two or more PDLs together along an arbitrary dimension (N-D
658 "append").
659 Sticks $x, $y, and all following arguments together along the specified
661 dimension. All other dimensions must be compatible in the broadcasting
662 sense.
663 Glue is permissive, in the sense that every PDL is treated as having an
665 infinite number of trivial dimensions of order 1 -- so "$x->glue(3,$y)"
666 works, even if $x and $y are only one dimensional.
667 If one of the PDLs has no elements, it is ignored. Likewise, if one of
669 them is actually the undefined value, it is treated as if it had no
670 elements.
671 If the first parameter is a defined perl scalar rather than a pdl, then
673 it is taken as a dimension along which to glue everything else, so you
674 can say "$cube = PDL::glue(3,@image_list);" if you like.
675 "glue" is implemented in pdl, using a combination of xchg and "append".
677 It should probably be updated (one day) to a pure PP function.
678 Similar functions include "append" (above), which appends only two
680 ndarrays along their first dimension, and cat, which can append more
681 than two ndarrays that all have the same sized dimensions.
682 srand
684 Signature: (a())
685 Seed random-number generator with a 64-bit int. Will generate seed data
687 for a number of threads equal to the return-value of "online_cpus" in
688 PDL::Core.
689 srand(); # uses current time
691 srand(5); # fixed number e.g. for testing
692 srand does not process bad values. It will set the bad-value flag of
694 all output ndarrays if the flag is set for any of the input ndarrays.
695 random
697 Signature: (a())
698 Constructor which returns ndarray of random numbers
700 $x = random([type], $nx, $ny, $nz,...);
702 $x = random $y;
703 etc (see zeroes).
705 This is the uniform distribution between 0 and 1 (assumedly excluding 1
707 itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
708 specify dimensions, types or give a template.
709 You can use the PDL function "srand" to seed the random generator. If
711 it has not been called yet, it will be with the current time.
712 random does not process bad values. It will set the bad-value flag of
714 all output ndarrays if the flag is set for any of the input ndarrays.
715 randsym
717 Signature: (a())
718 Constructor which returns ndarray of random numbers
720 $x = randsym([type], $nx, $ny, $nz,...);
722 $x = randsym $y;
723 etc (see zeroes).
725 This is the uniform distribution between 0 and 1 (excluding both 0 and
727 1, cf "random"). The arguments are the same as "zeroes" (q.v.) - i.e.
728 one can specify dimensions, types or give a template.
729 You can use the PDL function "srand" to seed the random generator. If
731 it has not been called yet, it will be with the current time.
732 randsym does not process bad values. It will set the bad-value flag of
734 all output ndarrays if the flag is set for any of the input ndarrays.
735 grandom
737 Constructor which returns ndarray of Gaussian random numbers
738 $x = grandom([type], $nx, $ny, $nz,...);
740 $x = grandom $y;
741 etc (see zeroes).
743 This is generated using the math library routine "ndtri".
745 Mean = 0, Stddev = 1
747 You can use the PDL function "srand" to seed the random generator. If
749 it has not been called yet, it will be with the current time.
750 vsearch
752 Signature: ( vals(); xs(n); [o] indx(); [\%options] )
753 Efficiently search for values in a sorted ndarray, returning indices.
755 $idx = vsearch( $vals, $x, [\%options] );
757 vsearch( $vals, $x, $idx, [\%options ] );
758 vsearch performs a binary search in the ordered ndarray $x, for the
760 values from $vals ndarray, returning indices into $x. What is a
761 "match", and the meaning of the returned indices, are determined by the
762 options.
763 The "mode" option indicates which method of searching to use, and may
765 be one of:
766 "sample"
768 invoke vsearch_sample, returning indices appropriate for sampling
769 within a distribution.
770 "insert_leftmost"
772 invoke vsearch_insert_leftmost, returning the left-most possible
773 insertion point which still leaves the ndarray sorted.
774 "insert_rightmost"
776 invoke vsearch_insert_rightmost, returning the right-most possible
777 insertion point which still leaves the ndarray sorted.
778 "match"
780 invoke vsearch_match, returning the index of a matching element,
781 else -(insertion point + 1)
782 "bin_inclusive"
784 invoke vsearch_bin_inclusive, returning an index appropriate for
785 binning on a grid where the left bin edges are inclusive of the
786 bin. See below for further explanation of the bin.
787 "bin_exclusive"
789 invoke vsearch_bin_exclusive, returning an index appropriate for
790 binning on a grid where the left bin edges are exclusive of the
791 bin. See below for further explanation of the bin.
792 The default value of "mode" is "sample".
794 use PDL;
796 my @modes = qw( sample insert_leftmost insert_rightmost match
798 bin_inclusive bin_exclusive );
799 # Generate a sequence of 3 zeros, 3 ones, ..., 3 fours.
801 my $x = zeroes(3,5)->yvals->flat;
802 for my $mode ( @modes ) {
804 # if the value is in $x
805 my $contained = 2;
806 my $idx_contained = vsearch( $contained, $x, { mode => $mode } );
807 my $x_contained = $x->copy;
808 $x_contained->slice( $idx_contained ) .= 9;
809 # if the value is not in $x
811 my $not_contained = 1.5;
812 my $idx_not_contained = vsearch( $not_contained, $x, { mode => $mode } );
813 my $x_not_contained = $x->copy;
814 $x_not_contained->slice( $idx_not_contained ) .= 9;
815 print sprintf("%-23s%30s\n", '$x', $x);
817 print sprintf("%-23s%30s\n", "$mode ($contained)", $x_contained);
818 print sprintf("%-23s%30s\n\n", "$mode ($not_contained)", $x_not_contained);
819 }
820 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
822 # sample (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
823 # sample (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
824 #
825 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
826 # insert_leftmost (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
827 # insert_leftmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
828 #
829 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
830 # insert_rightmost (2) [0 0 0 1 1 1 2 2 2 9 3 3 4 4 4]
831 # insert_rightmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
832 #
833 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
834 # match (2) [0 0 0 1 1 1 2 9 2 3 3 3 4 4 4]
835 # match (1.5) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
836 #
837 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
838 # bin_inclusive (2) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
839 # bin_inclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
840 #
841 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
842 # bin_exclusive (2) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
843 # bin_exclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
844 Also see vsearch_sample, vsearch_insert_leftmost,
846 vsearch_insert_rightmost, vsearch_match, vsearch_bin_inclusive, and
847 vsearch_bin_exclusive
848 vsearch_sample
850 Signature: (vals(); x(n); indx [o]idx())
851 Search for values in a sorted array, return index appropriate for
853 sampling from a distribution
854 $idx = vsearch_sample($vals, $x);
856 $x must be sorted, but may be in decreasing or increasing order.
858 vsearch_sample returns an index I for each value V of $vals appropriate
860 for sampling $vals
861 I has the following properties:
863 • if $x is sorted in increasing order
865 V <= x[0] : I = 0
867 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
868 x[-1] < V : I = $x->nelem -1
869 • if $x is sorted in decreasing order
871 V > x[0] : I = 0
873 x[0] >= V > x[-1] : I s.t. x[I] >= V > x[I+1]
874 x[-1] >= V : I = $x->nelem - 1
875 If all elements of $x are equal, I = $x->nelem - 1.
877 If $x contains duplicated elements, I is the index of the leftmost (by
879 position in array) duplicate if V matches.
880 This function is useful e.g. when you have a list of probabilities for
882 events and want to generate indices to events:
883 $x = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
885 $y = random 20;
886 $c = vsearch_sample($y, $x); # Now, $c will have the appropriate distr.
887 It is possible to use the cumusumover function to obtain cumulative
889 probabilities from absolute probabilities.
890 needs major (?) work to handles bad values
892 vsearch_insert_leftmost
894 Signature: (vals(); x(n); indx [o]idx())
895 Determine the insertion point for values in a sorted array, inserting
897 before duplicates.
898 $idx = vsearch_insert_leftmost($vals, $x);
900 $x must be sorted, but may be in decreasing or increasing order.
902 vsearch_insert_leftmost returns an index I for each value V of $vals
904 equal to the leftmost position (by index in array) within $x that V may
905 be inserted and still maintain the order in $x.
906 Insertion at index I involves shifting elements I and higher of $x to
908 the right by one and setting the now empty element at index I to V.
909 I has the following properties:
911 • if $x is sorted in increasing order
913 V <= x[0] : I = 0
915 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
916 x[-1] < V : I = $x->nelem
917 • if $x is sorted in decreasing order
919 V > x[0] : I = -1
921 x[0] >= V >= x[-1] : I s.t. x[I] >= V > x[I+1]
922 x[-1] >= V : I = $x->nelem -1
923 If all elements of $x are equal,
925 i = 0
927 If $x contains duplicated elements, I is the index of the leftmost (by
929 index in array) duplicate if V matches.
930 needs major (?) work to handles bad values
932 vsearch_insert_rightmost
934 Signature: (vals(); x(n); indx [o]idx())
935 Determine the insertion point for values in a sorted array, inserting
937 after duplicates.
938 $idx = vsearch_insert_rightmost($vals, $x);
940 $x must be sorted, but may be in decreasing or increasing order.
942 vsearch_insert_rightmost returns an index I for each value V of $vals
944 equal to the rightmost position (by index in array) within $x that V
945 may be inserted and still maintain the order in $x.
946 Insertion at index I involves shifting elements I and higher of $x to
948 the right by one and setting the now empty element at index I to V.
949 I has the following properties:
951 • if $x is sorted in increasing order
953 V < x[0] : I = 0
955 x[0] <= V < x[-1] : I s.t. x[I-1] <= V < x[I]
956 x[-1] <= V : I = $x->nelem
957 • if $x is sorted in decreasing order
959 V >= x[0] : I = -1
961 x[0] > V >= x[-1] : I s.t. x[I] >= V > x[I+1]
962 x[-1] > V : I = $x->nelem -1
963 If all elements of $x are equal,
965 i = $x->nelem - 1
967 If $x contains duplicated elements, I is the index of the leftmost (by
969 index in array) duplicate if V matches.
970 needs major (?) work to handles bad values
972 vsearch_match
974 Signature: (vals(); x(n); indx [o]idx())
975 Match values against a sorted array.
977 $idx = vsearch_match($vals, $x);
979 $x must be sorted, but may be in decreasing or increasing order.
981 vsearch_match returns an index I for each value V of $vals. If V
983 matches an element in $x, I is the index of that element, otherwise it
984 is -( insertion_point + 1 ), where insertion_point is an index in $x
985 where V may be inserted while maintaining the order in $x. If $x has
986 duplicated values, I may refer to any of them.
987 needs major (?) work to handles bad values
989 vsearch_bin_inclusive
991 Signature: (vals(); x(n); indx [o]idx())
992 Determine the index for values in a sorted array of bins, lower bound
994 inclusive.
995 $idx = vsearch_bin_inclusive($vals, $x);
997 $x must be sorted, but may be in decreasing or increasing order.
999 $x represents the edges of contiguous bins, with the first and last
1001 elements representing the outer edges of the outer bins, and the inner
1002 elements the shared bin edges.
1003 The lower bound of a bin is inclusive to the bin, its outer bound is
1005 exclusive to it. vsearch_bin_inclusive returns an index I for each
1006 value V of $vals
1007 I has the following properties:
1009 • if $x is sorted in increasing order
1011 V < x[0] : I = -1
1013 x[0] <= V < x[-1] : I s.t. x[I] <= V < x[I+1]
1014 x[-1] <= V : I = $x->nelem - 1
1015 • if $x is sorted in decreasing order
1017 V >= x[0] : I = 0
1019 x[0] > V >= x[-1] : I s.t. x[I+1] > V >= x[I]
1020 x[-1] > V : I = $x->nelem
1021 If all elements of $x are equal,
1023 i = $x->nelem - 1
1025 If $x contains duplicated elements, I is the index of the righmost (by
1027 index in array) duplicate if V matches.
1028 needs major (?) work to handles bad values
1030 vsearch_bin_exclusive
1032 Signature: (vals(); x(n); indx [o]idx())
1033 Determine the index for values in a sorted array of bins, lower bound
1035 exclusive.
1036 $idx = vsearch_bin_exclusive($vals, $x);
1038 $x must be sorted, but may be in decreasing or increasing order.
1040 $x represents the edges of contiguous bins, with the first and last
1042 elements representing the outer edges of the outer bins, and the inner
1043 elements the shared bin edges.
1044 The lower bound of a bin is exclusive to the bin, its upper bound is
1046 inclusive to it. vsearch_bin_exclusive returns an index I for each
1047 value V of $vals.
1048 I has the following properties:
1050 • if $x is sorted in increasing order
1052 V <= x[0] : I = -1
1054 x[0] < V <= x[-1] : I s.t. x[I] < V <= x[I+1]
1055 x[-1] < V : I = $x->nelem - 1
1056 • if $x is sorted in decreasing order
1058 V > x[0] : I = 0
1060 x[0] >= V > x[-1] : I s.t. x[I-1] >= V > x[I]
1061 x[-1] >= V : I = $x->nelem
1062 If all elements of $x are equal,
1064 i = $x->nelem - 1
1066 If $x contains duplicated elements, I is the index of the righmost (by
1068 index in array) duplicate if V matches.
1069 needs major (?) work to handles bad values
1071 interpolate
1073 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
1074 routine for 1D linear interpolation
1076 ( $yi, $err ) = interpolate($xi, $x, $y)
1078 Given a set of points "($x,$y)", use linear interpolation to find the
1080 values $yi at a set of points $xi.
1081 "interpolate" uses a binary search to find the suspects, er...,
1083 interpolation indices and therefore abscissas (ie $x) have to be
1084 strictly ordered (increasing or decreasing). For interpolation at lots
1085 of closely spaced abscissas an approach that uses the last index found
1086 as a start for the next search can be faster (compare Numerical Recipes
1087 "hunt" routine). Feel free to implement that on top of the binary
1088 search if you like. For out of bounds values it just does a linear
1089 extrapolation and sets the corresponding element of $err to 1, which is
1090 otherwise 0.
1091 See also "interpol", which uses the same routine, differing only in the
1093 handling of extrapolation - an error message is printed rather than
1094 returning an error ndarray.
1095 needs major (?) work to handles bad values
1097 interpol
1099 Signature: (xi(); x(n); y(n); [o] yi())
1100 routine for 1D linear interpolation
1102 $yi = interpol($xi, $x, $y)
1104 "interpol" uses the same search method as "interpolate", hence $x must
1106 be strictly ordered (either increasing or decreasing). The difference
1107 occurs in the handling of out-of-bounds values; here an error message
1108 is printed.
1109 interpND
1111 Interpolate values from an N-D ndarray, with switchable method
1112 $source = 10*xvals(10,10) + yvals(10,10);
1114 $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
1115 print $source->interpND( $index );
1116 InterpND acts like indexND, collapsing $index by lookup into $source;
1118 but it does interpolation rather than direct sampling. The
1119 interpolation method and boundary condition are switchable via an
1120 options hash.
1121 By default, linear or sample interpolation is used, with constant value
1123 outside the boundaries of the source pdl. No dataflow occurs, because
1124 in general the output is computed rather than indexed.
1125 All the interpolation methods treat the pixels as value-centered, so
1127 the "sample" method will return "$a->(0)" for coordinate values on the
1128 set [-0.5,0.5), and all methods will return "$a->(1)" for a coordinate
1129 value of exactly 1.
1130 Recognized options:
1132 method
1134 Values can be:
1135 • 0, s, sample, Sample (default for integer source types)
1137 The nearest value is taken. Pixels are regarded as centered on
1139 their respective integer coordinates (no offset from the linear
1140 case).
1141 • 1, l, linear, Linear (default for floating point source types)
1143 The values are N-linearly interpolated from an N-dimensional cube
1145 of size 2.
1146 • 3, c, cube, cubic, Cubic
1148 The values are interpolated using a local cubic fit to the data.
1150 The fit is constrained to match the original data and its
1151 derivative at the data points. The second derivative of the fit
1152 is not continuous at the data points. Multidimensional datasets
1153 are interpolated by the successive-collapse method.
1154 (Note that the constraint on the first derivative causes a small
1156 amount of ringing around sudden features such as step functions).
1157 • f, fft, fourier, Fourier
1159 The source is Fourier transformed, and the interpolated values
1161 are explicitly calculated from the coefficients. The boundary
1162 condition option is ignored -- periodic boundaries are imposed.
1163 If you pass in the option "fft", and it is a list (ARRAY) ref,
1165 then it is a stash for the magnitude and phase of the source FFT.
1166 If the list has two elements then they are taken as already
1167 computed; otherwise they are calculated and put in the stash.
1168 b, bound, boundary, Boundary
1170 This option is passed unmodified into indexND, which is used as the
1171 indexing engine for the interpolation. Some current allowed values
1172 are 'extend', 'periodic', 'truncate', and 'mirror' (default is
1173 'truncate').
1174 bad
1176 contains the fill value used for 'truncate' boundary. (default 0)
1177 fft
1179 An array ref whose associated list is used to stash the FFT of the
1180 source data, for the FFT method.
1181 one2nd
1183 Converts a one dimensional index ndarray to a set of ND coordinates
1184 @coords=one2nd($x, $indices)
1186 returns an array of ndarrays containing the ND indexes corresponding to
1188 the one dimensional list indices. The indices are assumed to correspond
1189 to array $x clumped using "clump(-1)". This routine is used in the old
1190 vector form of "whichND", but is useful on its own occasionally.
1191 Returned ndarrays have the indx datatype. $indices can have values
1193 larger than "$x->nelem" but negative values in $indices will not give
1194 the answer you expect.
1195 pdl> $x=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$x->clump(-1)
1197 pdl> $maxind=maximum_ind($c); p $maxind;
1198 6
1199 pdl> print one2nd($x, maximum_ind($c))
1200 0 1 1
1201 pdl> p $x->at(0,1,1)
1202 3
1203 which
1205 Signature: (mask(n); indx [o] inds(m))
1206 Returns indices of non-zero values from a 1-D PDL
1208 $i = which($mask);
1210 returns a pdl with indices for all those elements that are nonzero in
1212 the mask. Note that the returned indices will be 1D. If you feed in a
1213 multidimensional mask, it will be flattened before the indices are
1214 calculated. See also "whichND" for multidimensional masks.
1215 If you want to index into the original mask or a similar ndarray with
1217 output from "which", remember to flatten it before calling index:
1218 $data = random 5, 5;
1220 $idx = which $data > 0.5; # $idx is now 1D
1221 $bigsum = $data->flat->index($idx)->sum; # flatten before indexing
1222 Compare also "where" for similar functionality.
1224 SEE ALSO:
1226 "which_both" returns separately the indices of both zero and nonzero
1228 values in the mask.
1229 "where" returns associated values from a data PDL, rather than indices
1231 into the mask PDL.
1232 "whichND" returns N-D indices into a multidimensional PDL.
1234 pdl> $x = sequence(10); p $x
1236 [0 1 2 3 4 5 6 7 8 9]
1237 pdl> $indx = which($x>6); p $indx
1238 [7 8 9]
1239 which processes bad values. It will set the bad-value flag of all
1241 output ndarrays if the flag is set for any of the input ndarrays.
1242 which_both
1244 Signature: (mask(n); indx [o] inds(m); indx [o]notinds(q))
1245 Returns indices of zero and nonzero values in a mask PDL
1247 ($i, $c_i) = which_both($mask);
1249 This works just as "which", but the complement of $i will be in $c_i.
1251 pdl> $x = sequence(10); p $x
1253 [0 1 2 3 4 5 6 7 8 9]
1254 pdl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
1255 [5 6 7 8 9]
1256 [0 1 2 3 4]
1257 which_both processes bad values. It will set the bad-value flag of all
1259 output ndarrays if the flag is set for any of the input ndarrays.
1260 where
1262 Use a mask to select values from one or more data PDLs
1263 "where" accepts one or more data ndarrays and a mask ndarray. It
1265 returns a list of output ndarrays, corresponding to the input data
1266 ndarrays. Each output ndarray is a 1-dimensional list of values in its
1267 corresponding data ndarray. The values are drawn from locations where
1268 the mask is nonzero.
1269 The output PDLs are still connected to the original data PDLs, for the
1271 purpose of dataflow.
1272 "where" combines the functionality of "which" and index into a single
1274 operation.
1275 BUGS:
1277 While "where" works OK for most N-dimensional cases, it does not
1279 broadcast properly over (for example) the (N+1)th dimension in data
1280 that is compared to an N-dimensional mask. Use "whereND" for that.
1281 $i = $x->where($x+5 > 0); # $i contains those elements of $x
1283 # where mask ($x+5 > 0) is 1
1284 $i .= -5; # Set those elements (of $x) to -5. Together, these
1285 # commands clamp $x to a maximum of -5.
1286 It is also possible to use the same mask for several ndarrays with the
1288 same call:
1289 ($i,$j,$k) = where($x,$y,$z, $x+5>0);
1291 Note: $i is always 1-D, even if $x is >1-D.
1293 WARNING: The first argument (the values) and the second argument (the
1295 mask) currently have to have the exact same dimensions (or horrible
1296 things happen). You *cannot* broadcast over a smaller mask, for
1297 example.
1298 whereND
1300 "where" with support for ND masks and broadcasting
1301 "whereND" accepts one or more data ndarrays and a mask ndarray. It
1303 returns a list of output ndarrays, corresponding to the input data
1304 ndarrays. The values are drawn from locations where the mask is
1305 nonzero.
1306 "whereND" differs from "where" in that the mask dimensionality is
1308 preserved which allows for proper broadcasting of the selection
1309 operation over higher dimensions.
1310 As with "where" the output PDLs are still connected to the original
1312 data PDLs, for the purpose of dataflow.
1313 $sdata = whereND $data, $mask
1315 ($s1, $s2, ..., $sn) = whereND $d1, $d2, ..., $dn, $mask
1316 where
1318 $data is M dimensional
1320 $mask is N < M dimensional
1321 dims($data) 1..N == dims($mask) 1..N
1322 with broadcasting over N+1 to M dimensions
1323 $data = sequence(4,3,2); # example data array
1325 $mask4 = (random(4)>0.5); # example 1-D mask array, has $n4 true values
1326 $mask43 = (random(4,3)>0.5); # example 2-D mask array, has $n43 true values
1327 $sdat4 = whereND $data, $mask4; # $sdat4 is a [$n4,3,2] pdl
1328 $sdat43 = whereND $data, $mask43; # $sdat43 is a [$n43,2] pdl
1329 Just as with "where", you can use the returned value in an assignment.
1331 That means that both of these examples are valid:
1332 # Used to create a new slice stored in $sdat4:
1334 $sdat4 = $data->whereND($mask4);
1335 $sdat4 .= 0;
1336 # Used in lvalue context:
1337 $data->whereND($mask4) .= 0;
1338 SEE ALSO:
1340 "whichND" returns N-D indices into a multidimensional PDL, from a mask.
1342 whichND
1344 Return the coordinates of non-zero values in a mask.
1345 WhichND returns the N-dimensional coordinates of each nonzero value in
1347 a mask PDL with any number of dimensions. The returned values arrive
1348 as an array-of-vectors suitable for use in indexND or range.
1349 $coords = whichND($mask);
1351 returns a PDL containing the coordinates of the elements that are non-
1353 zero in $mask, suitable for use in "indexND" in PDL::Slices. The 0th
1354 dimension contains the full coordinate listing of each point; the 1st
1355 dimension lists all the points. For example, if $mask has rank 4 and
1356 100 matching elements, then $coords has dimension 4x100.
1357 If no such elements exist, then whichND returns a structured empty PDL:
1359 an Nx0 PDL that contains no values (but matches, broadcasting-wise,
1360 with the vectors that would be produced if such elements existed).
1361 DEPRECATED BEHAVIOR IN LIST CONTEXT:
1363 whichND once delivered different values in list context than in scalar
1365 context, for historical reasons. In list context, it returned the
1366 coordinates transposed, as a collection of 1-PDLs (one per dimension)
1367 in a list. This usage is deprecated in PDL 2.4.10, and will cause a
1368 warning to be issued every time it is encountered. To avoid the
1369 warning, you can set the global variable "$PDL::whichND" to 's' to get
1370 scalar behavior in all contexts, or to 'l' to get list behavior in list
1371 context.
1372 In later versions of PDL, the deprecated behavior will disappear.
1374 Deprecated list context whichND expressions can be replaced with:
1375 @list = $x->whichND->mv(0,-1)->dog;
1377 SEE ALSO:
1379 "which" finds coordinates of nonzero values in a 1-D mask.
1381 "where" extracts values from a data PDL that are associated with
1383 nonzero values in a mask PDL.
1384 "indexND" in PDL::Slices can be fed the coordinates to return the
1386 values.
1387 pdl> $s=sequence(10,10,3,4)
1389 pdl> ($x, $y, $z, $w)=whichND($s == 203); p $x, $y, $z, $w
1390 [3] [0] [2] [0]
1391 pdl> print $s->at(list(cat($x,$y,$z,$w)))
1392 203
1393 setops
1395 Implements simple set operations like union and intersection
1396 Usage: $set = setops($x, <OPERATOR>, $y);
1398 The operator can be "OR", "XOR" or "AND". This is then applied to $x
1400 viewed as a set and $y viewed as a set. Set theory says that a set may
1401 not have two or more identical elements, but setops takes care of this
1402 for you, so "$x=pdl(1,1,2)" is OK. The functioning is as follows:
1403 "OR"
1405 The resulting vector will contain the elements that are either in
1406 $x or in $y or both. This is the union in set operation terms
1407 "XOR"
1409 The resulting vector will contain the elements that are either in
1410 $x or $y, but not in both. This is
1411 Union($x, $y) - Intersection($x, $y)
1413 in set operation terms.
1415 "AND"
1417 The resulting vector will contain the intersection of $x and $y, so
1418 the elements that are in both $x and $y. Note that for convenience
1419 this operation is also aliased to "intersect".
1420 It should be emphasized that these routines are used when one or both
1422 of the sets $x, $y are hard to calculate or that you get from a
1423 separate subroutine.
1424 Finally IDL users might be familiar with Craig Markwardt's
1426 "cmset_op.pro" routine which has inspired this routine although it was
1427 written independently However the present routine has a few less
1428 options (but see the examples)
1429 You will very often use these functions on an index vector, so that is
1431 what we will show here. We will in fact something slightly silly. First
1432 we will find all squares that are also cubes below 10000.
1433 Create a sequence vector:
1435 pdl> $x = sequence(10000)
1437 Find all odd and even elements:
1439 pdl> ($even, $odd) = which_both( ($x % 2) == 0)
1441 Find all squares
1443 pdl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
1445 Find all cubes (being careful with roundoff error!)
1447 pdl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
1449 Then find all squares that are cubes:
1451 pdl> $both = setops($squares, 'AND', $cubes)
1453 And print these (assumes that "PDL::NiceSlice" is loaded!)
1455 pdl> p $x($both)
1457 [0 1 64 729 4096]
1458 Then find all numbers that are either cubes or squares, but not both:
1460 pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)
1462 pdl> p $cube_xor_square->nelem()
1464 112
1465 So there are a total of 112 of these!
1467 Finally find all odd squares:
1469 pdl> $odd_squares = setops($squares, 'AND', $odd)
1471 Another common occurrence is to want to get all objects that are in $x
1473 and in the complement of $y. But it is almost always best to create the
1474 complement explicitly since the universe that both are taken from is
1475 not known. Thus use "which_both" if possible to keep track of
1476 complements.
1477 If this is impossible the best approach is to make a temporary:
1479 This creates an index vector the size of the universe of the sets and
1481 set all elements in $y to 0
1482 pdl> $tmp = ones($n_universe); $tmp($y) .= 0;
1484 This then finds the complement of $y
1486 pdl> $C_b = which($tmp == 1);
1488 and this does the final selection:
1490 pdl> $set = setops($x, 'AND', $C_b)
1492 intersect
1494 Calculate the intersection of two ndarrays
1495 Usage: $set = intersect($x, $y);
1497 This routine is merely a simple interface to "setops". See that for
1499 more information
1500 Find all numbers less that 100 that are of the form 2*y and 3*x
1502 pdl> $x=sequence(100)
1504 pdl> $factor2 = which( ($x % 2) == 0)
1505 pdl> $factor3 = which( ($x % 3) == 0)
1506 pdl> $ii=intersect($factor2, $factor3)
1507 pdl> p $x($ii)
1508 [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
1509
1511 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
1512 Contributions by Christian Soeller (c.soeller@auckland.ac.nz), Karl
1513 Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
1514 (deforest@boulder.swri.edu) and Jarle Brinchmann (jarle@astro.up.pt)
1515 All rights reserved. There is no warranty. You are allowed to
1516 redistribute this software / documentation under certain conditions.
1517 For details, see the file COPYING in the PDL distribution. If this file
1518 is separated from the PDL distribution, the copyright notice should be
1519 included in the file.
1520 Updated for CPAN viewing compatibility by David Mertens.
1522perl v5.34.0 2022-02-28 Primitive(3)