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 threading 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 piddles if the flag is set for any of the input piddles.
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 Threading 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 thread over either of those dims.
63 EXAMPLES
65 Here are some simple ways to define vectors and matrices:
67 pdl> $r = pdl(1,2); # A row vector
69 pdl> $c = pdl([[3],[4]]); # A column vector
70 pdl> $c = pdl(3,4)->(*1); # A column vector, using NiceSlice
71 pdl> $m = pdl([[1,2],[3,4]]); # A 2x2 matrix
72 Now that we have a few objects prepared, here is how to matrix-multiply
74 them:
75 pdl> print $r x $m # row x matrix = row
77 [
78 [ 7 10]
79 ]
80 pdl> print $m x $r # matrix x row = ERROR
82 PDL: Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1
83 pdl> print $m x $c # matrix x column = column
85 [
86 [ 5]
87 [11]
88 ]
89 pdl> print $m x 2 # Trivial case: scalar mult.
91 [
92 [2 4]
93 [6 8]
94 ]
95 pdl> print $r x $c # row x column = scalar
97 [
98 [11]
99 ]
100 pdl> print $c x $r # column x row = matrix
102 [
103 [3 6]
104 [4 8]
105 ]
106 INTERNALS
108 The mechanics of the multiplication are carried out by the "matmult"
110 method.
111 matmult
113 Signature: (a(t,h); b(w,t); [o]c(w,h))
114 Matrix multiplication
116 Notionally, matrix multiplication $x x $y is equivalent to the
118 threading expression
119 $x->dummy(1)->inner($y->xchg(0,1)->dummy(2),$c);
121 but for large matrices that breaks CPU cache and is slow. Instead,
123 matmult calculates its result in 32x32x32 tiles, to keep the memory
124 footprint within cache as long as possible on most modern CPUs.
125 For usage, see "x", a description of the overloaded 'x' operator
127 matmult ignores the bad-value flag of the input piddles. It will set
129 the bad-value flag of all output piddles if the flag is set for any of
130 the input piddles.
131 innerwt
133 Signature: (a(n); b(n); c(n); [o]d())
134 Weighted (i.e. triple) inner product
136 d = sum_i a(i) b(i) c(i)
138 innerwt processes bad values. It will set the bad-value flag of all
140 output piddles if the flag is set for any of the input piddles.
141 inner2
143 Signature: (a(n); b(n,m); c(m); [o]d())
144 Inner product of two vectors and a matrix
146 d = sum_ij a(i) b(i,j) c(j)
148 Note that you should probably not thread over "a" and "c" since that
150 would be very wasteful. Instead, you should use a temporary for "b*c".
151 inner2 processes bad values. It will set the bad-value flag of all
153 output piddles if the flag is set for any of the input piddles.
154 inner2d
156 Signature: (a(n,m); b(n,m); [o]c())
157 Inner product over 2 dimensions.
159 Equivalent to
161 $c = inner($x->clump(2), $y->clump(2))
163 inner2d processes bad values. It will set the bad-value flag of all
165 output piddles if the flag is set for any of the input piddles.
166 inner2t
168 Signature: (a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k)))
169 Efficient Triple matrix product "a*b*c"
171 Efficiency comes from by using the temporary "tmp". This operation only
173 scales as "N**3" whereas threading using "inner2" would scale as
174 "N**4".
175 The reason for having this routine is that you do not need to have the
177 same thread-dimensions for "tmp" as for the other arguments, which in
178 case of large numbers of matrices makes this much more memory-
179 efficient.
180 It is hoped that things like this could be taken care of as a kind of
182 closures at some point.
183 inner2t processes bad values. It will set the bad-value flag of all
185 output piddles if the flag is set for any of the input piddles.
186 crossp
188 Signature: (a(tri=3); b(tri); [o] c(tri))
189 Cross product of two 3D vectors
191 After
193 $c = crossp $x, $y
195 the inner product "$c*$x" and "$c*$y" will be zero, i.e. $c is
197 orthogonal to $x and $y
198 crossp does not process bad values. It will set the bad-value flag of
200 all output piddles if the flag is set for any of the input piddles.
201 norm
203 Signature: (vec(n); [o] norm(n))
204 Normalises a vector to unit Euclidean length
206 norm processes bad values. It will set the bad-value flag of all
208 output piddles if the flag is set for any of the input piddles.
209 indadd
211 Signature: (a(); indx ind(); [o] sum(m))
212 Threaded Index Add: Add "a" to the "ind" element of "sum", i.e:
214 sum(ind) += a
216 Simple Example:
218 $x = 2;
220 $ind = 3;
221 $sum = zeroes(10);
222 indadd($x,$ind, $sum);
223 print $sum
224 #Result: ( 2 added to element 3 of $sum)
225 # [0 0 0 2 0 0 0 0 0 0]
226 Threaded Example:
228 $x = pdl( 1,2,3);
230 $ind = pdl( 1,4,6);
231 $sum = zeroes(10);
232 indadd($x,$ind, $sum);
233 print $sum."\n";
234 #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
235 # [0 1 0 0 2 0 3 0 0 0]
236 The routine barfs if any of the indices are bad.
238 conv1d
240 Signature: (a(m); kern(p); [o]b(m); int reflect)
241 1D convolution along first dimension
243 The m-th element of the discrete convolution of an input piddle $a of
245 size $M, and a kernel piddle $kern of size $P, is calculated as
246 n = ($P-1)/2
248 ====
249 \
250 ($a conv1d $kern)[m] = > $a_ext[m - n] * $kern[n]
251 /
252 ====
253 n = -($P-1)/2
254 where $a_ext is either the periodic (or reflected) extension of $a so
256 it is equal to $a on " 0..$M-1 " and equal to the corresponding
257 periodic/reflected image of $a outside that range.
258 $con = conv1d sequence(10), pdl(-1,0,1);
260 $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
262 By default, periodic boundary conditions are assumed (i.e. wrap
264 around). Alternatively, you can request reflective boundary conditions
265 using the "Boundary" option:
266 {Boundary => 'reflect'} # case in 'reflect' doesn't matter
268 The convolution is performed along the first dimension. To apply it
270 across another dimension use the slicing routines, e.g.
271 $y = $x->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
273 This function is useful for threaded filtering of 1D signals.
275 Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv
277 WARNING: "conv1d" processes bad values in its inputs as the numeric
279 value of "$pdl->badvalue" so it is not recommended for processing pdls
280 with bad values in them unless special care is taken.
281 conv1d ignores the bad-value flag of the input piddles. It will set
283 the bad-value flag of all output piddles if the flag is set for any of
284 the input piddles.
285 in
287 Signature: (a(); b(n); [o] c())
288 test if a is in the set of values b
290 $goodmsk = $labels->in($goodlabels);
292 print pdl(3,1,4,6,2)->in(pdl(2,3,3));
293 [1 0 0 0 1]
294 "in" is akin to the is an element of of set theory. In principle, PDL
296 threading could be used to achieve its functionality by using a
297 construct like
298 $msk = ($labels->dummy(0) == $goodlabels)->orover;
300 However, "in" doesn't create a (potentially large) intermediate and is
302 generally faster.
303 in does not process bad values. It will set the bad-value flag of all
305 output piddles if the flag is set for any of the input piddles.
306 uniq
308 return all unique elements of a piddle
309 The unique elements are returned in ascending order.
311 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniq
313 [-1 0 2 4 6] # 0 is returned 2nd (sorted order)
314 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniq
316 [-1 2 4 6 nan] # NaN value is returned at end
317 Note: The returned pdl is 1D; any structure of the input piddle is
319 lost. "NaN" values are never compare equal to any other values, even
320 themselves. As a result, they are always unique. "uniq" returns the
321 NaN values at the end of the result piddle. This follows the Matlab
322 usage.
323 See "uniqind" if you need the indices of the unique elements rather
325 than the values.
326 Bad values are not considered unique by uniq and are ignored.
328 $x=sequence(10);
330 $x=$x->setbadif($x%3);
331 print $x->uniq;
332 [0 3 6 9]
333 uniqind
335 Return the indices of all unique elements of a piddle The order is in
336 the order of the values to be consistent with uniq. "NaN" values never
337 compare equal with any other value and so are always unique. This
338 follows the Matlab usage.
339 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniqind
341 [5 4 1 3 6] # the 0 at index 4 is returned 2nd, but...
342 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniqind
344 [5 1 3 6 4] # ...the NaN at index 4 is returned at end
345 Note: The returned pdl is 1D; any structure of the input piddle is
347 lost.
348 See "uniq" if you want the unique values instead of the indices.
350 Bad values are not considered unique by uniqind and are ignored.
352 uniqvec
354 Return all unique vectors out of a collection
355 NOTE: If any vectors in the input piddle have NaN values
357 they are returned at the end of the non-NaN ones. This is
358 because, by definition, NaN values never compare equal with
359 any other value.
360 NOTE: The current implementation does not sort the vectors
362 containing NaN values.
363 The unique vectors are returned in lexicographically sorted ascending
365 order. The 0th dimension of the input PDL is treated as a dimensional
366 index within each vector, and the 1st and any higher dimensions are
367 taken to run across vectors. The return value is always 2D; any
368 structure of the input PDL (beyond using the 0th dimension for vector
369 index) is lost.
370 See also "uniq" for a unique list of scalars; and qsortvec for sorting
372 a list of vectors lexicographcally.
373 If a vector contains all bad values, it is ignored as in "uniq". If
375 some of the values are good, it is treated as a normal vector. For
376 example, [1 2 BAD] and [BAD 2 3] could be returned, but [BAD BAD BAD]
377 could not. Vectors containing BAD values will be returned after any
378 non-NaN and non-BAD containing vectors, followed by the NaN vectors.
379 hclip
381 Signature: (a(); b(); [o] c())
382 clip (threshold) $a by $b ($b is upper bound)
384 hclip processes bad values. It will set the bad-value flag of all
386 output piddles if the flag is set for any of the input piddles.
387 lclip
389 Signature: (a(); b(); [o] c())
390 clip (threshold) $a by $b ($b is lower bound)
392 lclip processes bad values. It will set the bad-value flag of all
394 output piddles if the flag is set for any of the input piddles.
395 clip
397 Clip (threshold) a piddle by (optional) upper or lower bounds.
398 $y = $x->clip(0,3);
400 $c = $x->clip(undef, $x);
401 clip handles bad values since it is just a wrapper around "hclip" and
403 "lclip".
404 clip
406 Signature: (a(); l(); h(); [o] c())
407 info not available
409 clip processes bad values. It will set the bad-value flag of all
411 output piddles if the flag is set for any of the input piddles.
412 wtstat
414 Signature: (a(n); wt(n); avg(); [o]b(); int deg)
415 Weighted statistical moment of given degree
417 This calculates a weighted statistic over the vector "a". The formula
419 is
420 b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
422 Bad values are ignored in any calculation; $b will only have its bad
424 flag set if the output contains any bad data.
425 statsover
427 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())
428 Calculate useful statistics over a dimension of a piddle
430 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($piddle, $weights);
432 This utility function calculates various useful quantities of a piddle.
434 These are:
435 • the mean:
437 MEAN = sum (x)/ N
439 with "N" being the number of elements in x
441 • the population RMS deviation from the mean:
443 PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
445 The population deviation is the best-estimate of the deviation of
447 the population from which a sample is drawn.
448 • the median
450 The median is the 50th percentile data value. Median is found by
452 medover, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
453 • the minimum
455 • the maximum
457 • the average absolute deviation:
459 AADEV = sum( abs(x-mean(x)) )/N
461 • RMS deviation from the mean:
463 RMS = sqrt(sum( (x-mean(x))^2 )/N)
465 (also known as the root-mean-square deviation, or the square root of
467 the variance)
468 This operator is a projection operator so the calculation will take
470 place over the final dimension. Thus if the input is N-dimensional each
471 returned value will be N-1 dimensional, to calculate the statistics for
472 the entire piddle either use "clump(-1)" directly on the piddle or call
473 "stats".
474 Bad values are simply ignored in the calculation, effectively reducing
476 the sample size. If all data are bad then the output data are marked
477 bad.
478 stats
480 Calculates useful statistics on a piddle
481 ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($piddle,[$weights]);
483 This utility calculates all the most useful quantities in one call. It
485 works the same way as "statsover", except that the quantities are
486 calculated considering the entire input PDL as a single sample, rather
487 than as a collection of rows. See "statsover" for definitions of the
488 returned quantities.
489 Bad values are handled; if all input values are bad, then all of the
491 output values are flagged bad.
492 histogram
494 Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
495 Calculates a histogram for given stepsize and minimum.
497 $h = histogram($data, $step, $min, $numbins);
499 $hist = zeroes $numbins; # Put histogram in existing piddle.
500 histogram($data, $hist, $step, $min, $numbins);
501 The histogram will contain $numbins bins starting from $min, each $step
503 wide. The value in each bin is the number of values in $data that lie
504 within the bin limits.
505 Data below the lower limit is put in the first bin, and data above the
507 upper limit is put in the last bin.
508 The output is reset in a different threadloop so that you can take a
510 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
511 For a higher-level interface, see hist.
513 pdl> p histogram(pdl(1,1,2),1,0,3)
515 [0 2 1]
516 histogram processes bad values. It will set the bad-value flag of all
518 output piddles if the flag is set for any of the input piddles.
519 whistogram
521 Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
522 Calculates a histogram from weighted data for given stepsize and
524 minimum.
525 $h = whistogram($data, $weights, $step, $min, $numbins);
527 $hist = zeroes $numbins; # Put histogram in existing piddle.
528 whistogram($data, $weights, $hist, $step, $min, $numbins);
529 The histogram will contain $numbins bins starting from $min, each $step
531 wide. The value in each bin is the sum of the values in $weights that
532 correspond to values in $data that lie within the bin limits.
533 Data below the lower limit is put in the first bin, and data above the
535 upper limit is put in the last bin.
536 The output is reset in a different threadloop so that you can take a
538 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
539 pdl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
541 [0 0.2 0.5 0]
542 whistogram processes bad values. It will set the bad-value flag of all
544 output piddles if the flag is set for any of the input piddles.
545 histogram2d
547 Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
548 double stepb; double minb; int mbsize => mb;)
549 Calculates a 2d histogram.
551 $h = histogram2d($datax, $datay, $stepx, $minx,
553 $nbinx, $stepy, $miny, $nbiny);
554 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
555 histogram2d($datax, $datay, $hist, $stepx, $minx,
556 $nbinx, $stepy, $miny, $nbiny);
557 The histogram will contain $nbinx x $nbiny bins, with the lower limits
559 of the first one at "($minx, $miny)", and with bin size "($stepx,
560 $stepy)". The value in each bin is the number of values in $datax and
561 $datay that lie within the bin limits.
562 Data below the lower limit is put in the first bin, and data above the
564 upper limit is put in the last bin.
565 pdl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
567 [
568 [0 0 0]
569 [0 2 2]
570 [0 1 0]
571 ]
572 histogram2d processes bad values. It will set the bad-value flag of
574 all output piddles if the flag is set for any of the input piddles.
575 whistogram2d
577 Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
578 double stepb; double minb; int mbsize => mb;)
579 Calculates a 2d histogram from weighted data.
581 $h = whistogram2d($datax, $datay, $weights,
583 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
584 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
585 whistogram2d($datax, $datay, $weights, $hist,
586 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
587 The histogram will contain $nbinx x $nbiny bins, with the lower limits
589 of the first one at "($minx, $miny)", and with bin size "($stepx,
590 $stepy)". The value in each bin is the sum of the values in $weights
591 that correspond to values in $datax and $datay that lie within the bin
592 limits.
593 Data below the lower limit is put in the first bin, and data above the
595 upper limit is put in the last bin.
596 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)
598 [
599 [ 0 0 0]
600 [ 0 0.5 0.9]
601 [ 0 0.1 0]
602 ]
603 whistogram2d processes bad values. It will set the bad-value flag of
605 all output piddles if the flag is set for any of the input piddles.
606 fibonacci
608 Signature: ([o]x(n))
609 Constructor - a vector with Fibonacci's sequence
611 fibonacci does not process bad values. It will set the bad-value flag
613 of all output piddles if the flag is set for any of the input piddles.
614 append
616 Signature: (a(n); b(m); [o] c(mn))
617 append two piddles by concatenating along their first dimensions
619 $x = ones(2,4,7);
621 $y = sequence 5;
622 $c = $x->append($y); # size of $c is now (7,4,7) (a jumbo-piddle ;)
623 "append" appends two piddles along their first dimensions. The rest of
625 the dimensions must be compatible in the threading sense. The resulting
626 size of the first dimension is the sum of the sizes of the first
627 dimensions of the two argument piddles - i.e. "n + m".
628 Similar functions include "glue" (below), which can append more than
630 two piddles along an arbitrary dimension, and cat, which can append
631 more than two piddles that all have the same sized dimensions.
632 append does not process bad values. It will set the bad-value flag of
634 all output piddles if the flag is set for any of the input piddles.
635 glue
637 $c = $x->glue(<dim>,$y,...)
638 Glue two or more PDLs together along an arbitrary dimension (N-D
640 "append").
641 Sticks $x, $y, and all following arguments together along the specified
643 dimension. All other dimensions must be compatible in the threading
644 sense.
645 Glue is permissive, in the sense that every PDL is treated as having an
647 infinite number of trivial dimensions of order 1 -- so "$x->glue(3,$y)"
648 works, even if $x and $y are only one dimensional.
649 If one of the PDLs has no elements, it is ignored. Likewise, if one of
651 them is actually the undefined value, it is treated as if it had no
652 elements.
653 If the first parameter is a defined perl scalar rather than a pdl, then
655 it is taken as a dimension along which to glue everything else, so you
656 can say "$cube = PDL::glue(3,@image_list);" if you like.
657 "glue" is implemented in pdl, using a combination of xchg and "append".
659 It should probably be updated (one day) to a pure PP function.
660 Similar functions include "append" (above), which appends only two
662 piddles along their first dimension, and cat, which can append more
663 than two piddles that all have the same sized dimensions.
664 axisvalues
666 Signature: ([o,nc]a(n))
667 Internal routine
669 "axisvalues" is the internal primitive that implements axisvals and
671 alters its argument.
672 axisvalues does not process bad values. It will set the bad-value flag
674 of all output piddles if the flag is set for any of the input piddles.
675 random
677 Constructor which returns piddle of random numbers
678 $x = random([type], $nx, $ny, $nz,...);
680 $x = random $y;
681 etc (see zeroes).
683 This is the uniform distribution between 0 and 1 (assumedly excluding 1
685 itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
686 specify dimensions, types or give a template.
687 You can use the perl function srand to seed the random generator. For
689 further details consult Perl's srand documentation.
690 randsym
692 Constructor which returns piddle of random numbers
693 $x = randsym([type], $nx, $ny, $nz,...);
695 $x = randsym $y;
696 etc (see zeroes).
698 This is the uniform distribution between 0 and 1 (excluding both 0 and
700 1, cf "random"). The arguments are the same as "zeroes" (q.v.) - i.e.
701 one can specify dimensions, types or give a template.
702 You can use the perl function srand to seed the random generator. For
704 further details consult Perl's srand documentation.
705 grandom
707 Constructor which returns piddle of Gaussian random numbers
708 $x = grandom([type], $nx, $ny, $nz,...);
710 $x = grandom $y;
711 etc (see zeroes).
713 This is generated using the math library routine "ndtri".
715 Mean = 0, Stddev = 1
717 You can use the perl function srand to seed the random generator. For
719 further details consult Perl's srand documentation.
720 vsearch
722 Signature: ( vals(); xs(n); [o] indx(); [\%options] )
723 Efficiently search for values in a sorted piddle, returning indices.
725 $idx = vsearch( $vals, $x, [\%options] );
727 vsearch( $vals, $x, $idx, [\%options ] );
728 vsearch performs a binary search in the ordered piddle $x, for the
730 values from $vals piddle, returning indices into $x. What is a
731 "match", and the meaning of the returned indices, are determined by the
732 options.
733 The "mode" option indicates which method of searching to use, and may
735 be one of:
736 "sample"
738 invoke vsearch_sample, returning indices appropriate for sampling
739 within a distribution.
740 "insert_leftmost"
742 invoke vsearch_insert_leftmost, returning the left-most possible
743 insertion point which still leaves the piddle sorted.
744 "insert_rightmost"
746 invoke vsearch_insert_rightmost, returning the right-most possible
747 insertion point which still leaves the piddle sorted.
748 "match"
750 invoke vsearch_match, returning the index of a matching element,
751 else -(insertion point + 1)
752 "bin_inclusive"
754 invoke vsearch_bin_inclusive, returning an index appropriate for
755 binning on a grid where the left bin edges are inclusive of the
756 bin. See below for further explanation of the bin.
757 "bin_exclusive"
759 invoke vsearch_bin_exclusive, returning an index appropriate for
760 binning on a grid where the left bin edges are exclusive of the
761 bin. See below for further explanation of the bin.
762 The default value of "mode" is "sample".
764 use PDL;
766 my @modes = qw( sample insert_leftmost insert_rightmost match
768 bin_inclusive bin_exclusive );
769 # Generate a sequence of 3 zeros, 3 ones, ..., 3 fours.
771 my $x = zeroes(3,5)->yvals->flat;
772 for my $mode ( @modes ) {
774 # if the value is in $x
775 my $contained = 2;
776 my $idx_contained = vsearch( $contained, $x, { mode => $mode } );
777 my $x_contained = $x->copy;
778 $x_contained->slice( $idx_contained ) .= 9;
779 # if the value is not in $x
781 my $not_contained = 1.5;
782 my $idx_not_contained = vsearch( $not_contained, $x, { mode => $mode } );
783 my $x_not_contained = $x->copy;
784 $x_not_contained->slice( $idx_not_contained ) .= 9;
785 print sprintf("%-23s%30s\n", '$x', $x);
787 print sprintf("%-23s%30s\n", "$mode ($contained)", $x_contained);
788 print sprintf("%-23s%30s\n\n", "$mode ($not_contained)", $x_not_contained);
789 }
790 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
792 # sample (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
793 # sample (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
794 #
795 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
796 # insert_leftmost (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
797 # insert_leftmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
798 #
799 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
800 # insert_rightmost (2) [0 0 0 1 1 1 2 2 2 9 3 3 4 4 4]
801 # insert_rightmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
802 #
803 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
804 # match (2) [0 0 0 1 1 1 2 9 2 3 3 3 4 4 4]
805 # match (1.5) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
806 #
807 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
808 # bin_inclusive (2) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
809 # bin_inclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
810 #
811 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
812 # bin_exclusive (2) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
813 # bin_exclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
814 Also see vsearch_sample, vsearch_insert_leftmost,
816 vsearch_insert_rightmost, vsearch_match, vsearch_bin_inclusive, and
817 vsearch_bin_exclusive
818 vsearch_sample
820 Signature: (vals(); x(n); indx [o]idx())
821 Search for values in a sorted array, return index appropriate for
823 sampling from a distribution
824 $idx = vsearch_sample($vals, $x);
826 $x must be sorted, but may be in decreasing or increasing order.
828 vsearch_sample returns an index I for each value V of $vals appropriate
830 for sampling $vals
831 I has the following properties:
833 • if $x is sorted in increasing order
835 V <= x[0] : I = 0
837 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
838 x[-1] < V : I = $x->nelem -1
839 • if $x is sorted in decreasing order
841 V > x[0] : I = 0
843 x[0] >= V > x[-1] : I s.t. x[I] >= V > x[I+1]
844 x[-1] >= V : I = $x->nelem - 1
845 If all elements of $x are equal, I = $x->nelem - 1.
847 If $x contains duplicated elements, I is the index of the leftmost (by
849 position in array) duplicate if V matches.
850 This function is useful e.g. when you have a list of probabilities for
852 events and want to generate indices to events:
853 $x = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
855 $y = random 20;
856 $c = vsearch_sample($y, $x); # Now, $c will have the appropriate distr.
857 It is possible to use the cumusumover function to obtain cumulative
859 probabilities from absolute probabilities.
860 needs major (?) work to handles bad values
862 vsearch_insert_leftmost
864 Signature: (vals(); x(n); indx [o]idx())
865 Determine the insertion point for values in a sorted array, inserting
867 before duplicates.
868 $idx = vsearch_insert_leftmost($vals, $x);
870 $x must be sorted, but may be in decreasing or increasing order.
872 vsearch_insert_leftmost returns an index I for each value V of $vals
874 equal to the leftmost position (by index in array) within $x that V may
875 be inserted and still maintain the order in $x.
876 Insertion at index I involves shifting elements I and higher of $x to
878 the right by one and setting the now empty element at index I to V.
879 I has the following properties:
881 • if $x is sorted in increasing order
883 V <= x[0] : I = 0
885 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
886 x[-1] < V : I = $x->nelem
887 • if $x is sorted in decreasing order
889 V > x[0] : I = -1
891 x[0] >= V >= x[-1] : I s.t. x[I] >= V > x[I+1]
892 x[-1] >= V : I = $x->nelem -1
893 If all elements of $x are equal,
895 i = 0
897 If $x contains duplicated elements, I is the index of the leftmost (by
899 index in array) duplicate if V matches.
900 needs major (?) work to handles bad values
902 vsearch_insert_rightmost
904 Signature: (vals(); x(n); indx [o]idx())
905 Determine the insertion point for values in a sorted array, inserting
907 after duplicates.
908 $idx = vsearch_insert_rightmost($vals, $x);
910 $x must be sorted, but may be in decreasing or increasing order.
912 vsearch_insert_rightmost returns an index I for each value V of $vals
914 equal to the rightmost position (by index in array) within $x that V
915 may be inserted and still maintain the order in $x.
916 Insertion at index I involves shifting elements I and higher of $x to
918 the right by one and setting the now empty element at index I to V.
919 I has the following properties:
921 • if $x is sorted in increasing order
923 V < x[0] : I = 0
925 x[0] <= V < x[-1] : I s.t. x[I-1] <= V < x[I]
926 x[-1] <= V : I = $x->nelem
927 • if $x is sorted in decreasing order
929 V >= x[0] : I = -1
931 x[0] > V >= x[-1] : I s.t. x[I] >= V > x[I+1]
932 x[-1] > V : I = $x->nelem -1
933 If all elements of $x are equal,
935 i = $x->nelem - 1
937 If $x contains duplicated elements, I is the index of the leftmost (by
939 index in array) duplicate if V matches.
940 needs major (?) work to handles bad values
942 vsearch_match
944 Signature: (vals(); x(n); indx [o]idx())
945 Match values against a sorted array.
947 $idx = vsearch_match($vals, $x);
949 $x must be sorted, but may be in decreasing or increasing order.
951 vsearch_match returns an index I for each value V of $vals. If V
953 matches an element in $x, I is the index of that element, otherwise it
954 is -( insertion_point + 1 ), where insertion_point is an index in $x
955 where V may be inserted while maintaining the order in $x. If $x has
956 duplicated values, I may refer to any of them.
957 needs major (?) work to handles bad values
959 vsearch_bin_inclusive
961 Signature: (vals(); x(n); indx [o]idx())
962 Determine the index for values in a sorted array of bins, lower bound
964 inclusive.
965 $idx = vsearch_bin_inclusive($vals, $x);
967 $x must be sorted, but may be in decreasing or increasing order.
969 $x represents the edges of contiguous bins, with the first and last
971 elements representing the outer edges of the outer bins, and the inner
972 elements the shared bin edges.
973 The lower bound of a bin is inclusive to the bin, its outer bound is
975 exclusive to it. vsearch_bin_inclusive returns an index I for each
976 value V of $vals
977 I has the following properties:
979 • if $x is sorted in increasing order
981 V < x[0] : I = -1
983 x[0] <= V < x[-1] : I s.t. x[I] <= V < x[I+1]
984 x[-1] <= V : I = $x->nelem - 1
985 • if $x is sorted in decreasing order
987 V >= x[0] : I = 0
989 x[0] > V >= x[-1] : I s.t. x[I+1] > V >= x[I]
990 x[-1] > V : I = $x->nelem
991 If all elements of $x are equal,
993 i = $x->nelem - 1
995 If $x contains duplicated elements, I is the index of the righmost (by
997 index in array) duplicate if V matches.
998 needs major (?) work to handles bad values
1000 vsearch_bin_exclusive
1002 Signature: (vals(); x(n); indx [o]idx())
1003 Determine the index for values in a sorted array of bins, lower bound
1005 exclusive.
1006 $idx = vsearch_bin_exclusive($vals, $x);
1008 $x must be sorted, but may be in decreasing or increasing order.
1010 $x represents the edges of contiguous bins, with the first and last
1012 elements representing the outer edges of the outer bins, and the inner
1013 elements the shared bin edges.
1014 The lower bound of a bin is exclusive to the bin, its upper bound is
1016 inclusive to it. vsearch_bin_exclusive returns an index I for each
1017 value V of $vals.
1018 I has the following properties:
1020 • if $x is sorted in increasing order
1022 V <= x[0] : I = -1
1024 x[0] < V <= x[-1] : I s.t. x[I] < V <= x[I+1]
1025 x[-1] < V : I = $x->nelem - 1
1026 • if $x is sorted in decreasing order
1028 V > x[0] : I = 0
1030 x[0] >= V > x[-1] : I s.t. x[I-1] >= V > x[I]
1031 x[-1] >= V : I = $x->nelem
1032 If all elements of $x are equal,
1034 i = $x->nelem - 1
1036 If $x contains duplicated elements, I is the index of the righmost (by
1038 index in array) duplicate if V matches.
1039 needs major (?) work to handles bad values
1041 interpolate
1043 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
1044 routine for 1D linear interpolation
1046 ( $yi, $err ) = interpolate($xi, $x, $y)
1048 Given a set of points "($x,$y)", use linear interpolation to find the
1050 values $yi at a set of points $xi.
1051 "interpolate" uses a binary search to find the suspects, er...,
1053 interpolation indices and therefore abscissas (ie $x) have to be
1054 strictly ordered (increasing or decreasing). For interpolation at lots
1055 of closely spaced abscissas an approach that uses the last index found
1056 as a start for the next search can be faster (compare Numerical Recipes
1057 "hunt" routine). Feel free to implement that on top of the binary
1058 search if you like. For out of bounds values it just does a linear
1059 extrapolation and sets the corresponding element of $err to 1, which is
1060 otherwise 0.
1061 See also "interpol", which uses the same routine, differing only in the
1063 handling of extrapolation - an error message is printed rather than
1064 returning an error piddle.
1065 needs major (?) work to handles bad values
1067 interpol
1069 Signature: (xi(); x(n); y(n); [o] yi())
1070 routine for 1D linear interpolation
1072 $yi = interpol($xi, $x, $y)
1074 "interpol" uses the same search method as "interpolate", hence $x must
1076 be strictly ordered (either increasing or decreasing). The difference
1077 occurs in the handling of out-of-bounds values; here an error message
1078 is printed.
1079 interpND
1081 Interpolate values from an N-D piddle, with switchable method
1082 $source = 10*xvals(10,10) + yvals(10,10);
1084 $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
1085 print $source->interpND( $index );
1086 InterpND acts like indexND, collapsing $index by lookup into $source;
1088 but it does interpolation rather than direct sampling. The
1089 interpolation method and boundary condition are switchable via an
1090 options hash.
1091 By default, linear or sample interpolation is used, with constant value
1093 outside the boundaries of the source pdl. No dataflow occurs, because
1094 in general the output is computed rather than indexed.
1095 All the interpolation methods treat the pixels as value-centered, so
1097 the "sample" method will return "$a->(0)" for coordinate values on the
1098 set [-0.5,0.5), and all methods will return "$a->(1)" for a coordinate
1099 value of exactly 1.
1100 Recognized options:
1102 method
1104 Values can be:
1105 • 0, s, sample, Sample (default for integer source types)
1107 The nearest value is taken. Pixels are regarded as centered on
1109 their respective integer coordinates (no offset from the linear
1110 case).
1111 • 1, l, linear, Linear (default for floating point source types)
1113 The values are N-linearly interpolated from an N-dimensional cube
1115 of size 2.
1116 • 3, c, cube, cubic, Cubic
1118 The values are interpolated using a local cubic fit to the data.
1120 The fit is constrained to match the original data and its
1121 derivative at the data points. The second derivative of the fit
1122 is not continuous at the data points. Multidimensional datasets
1123 are interpolated by the successive-collapse method.
1124 (Note that the constraint on the first derivative causes a small
1126 amount of ringing around sudden features such as step functions).
1127 • f, fft, fourier, Fourier
1129 The source is Fourier transformed, and the interpolated values
1131 are explicitly calculated from the coefficients. The boundary
1132 condition option is ignored -- periodic boundaries are imposed.
1133 If you pass in the option "fft", and it is a list (ARRAY) ref,
1135 then it is a stash for the magnitude and phase of the source FFT.
1136 If the list has two elements then they are taken as already
1137 computed; otherwise they are calculated and put in the stash.
1138 b, bound, boundary, Boundary
1140 This option is passed unmodified into indexND, which is used as the
1141 indexing engine for the interpolation. Some current allowed values
1142 are 'extend', 'periodic', 'truncate', and 'mirror' (default is
1143 'truncate').
1144 bad
1146 contains the fill value used for 'truncate' boundary. (default 0)
1147 fft
1149 An array ref whose associated list is used to stash the FFT of the
1150 source data, for the FFT method.
1151 one2nd
1153 Converts a one dimensional index piddle to a set of ND coordinates
1154 @coords=one2nd($x, $indices)
1156 returns an array of piddles containing the ND indexes corresponding to
1158 the one dimensional list indices. The indices are assumed to correspond
1159 to array $x clumped using "clump(-1)". This routine is used in the old
1160 vector form of "whichND", but is useful on its own occasionally.
1161 Returned piddles have the indx datatype. $indices can have values
1163 larger than "$x->nelem" but negative values in $indices will not give
1164 the answer you expect.
1165 pdl> $x=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$x->clump(-1)
1167 pdl> $maxind=maximum_ind($c); p $maxind;
1168 6
1169 pdl> print one2nd($x, maximum_ind($c))
1170 0 1 1
1171 pdl> p $x->at(0,1,1)
1172 3
1173 which
1175 Signature: (mask(n); indx [o] inds(m))
1176 Returns indices of non-zero values from a 1-D PDL
1178 $i = which($mask);
1180 returns a pdl with indices for all those elements that are nonzero in
1182 the mask. Note that the returned indices will be 1D. If you feed in a
1183 multidimensional mask, it will be flattened before the indices are
1184 calculated. See also "whichND" for multidimensional masks.
1185 If you want to index into the original mask or a similar piddle with
1187 output from "which", remember to flatten it before calling index:
1188 $data = random 5, 5;
1190 $idx = which $data > 0.5; # $idx is now 1D
1191 $bigsum = $data->flat->index($idx)->sum; # flatten before indexing
1192 Compare also "where" for similar functionality.
1194 SEE ALSO:
1196 "which_both" returns separately the indices of both zero and nonzero
1198 values in the mask.
1199 "where" returns associated values from a data PDL, rather than indices
1201 into the mask PDL.
1202 "whichND" returns N-D indices into a multidimensional PDL.
1204 pdl> $x = sequence(10); p $x
1206 [0 1 2 3 4 5 6 7 8 9]
1207 pdl> $indx = which($x>6); p $indx
1208 [7 8 9]
1209 which processes bad values. It will set the bad-value flag of all
1211 output piddles if the flag is set for any of the input piddles.
1212 which_both
1214 Signature: (mask(n); indx [o] inds(m); indx [o]notinds(q))
1215 Returns indices of zero and nonzero values in a mask PDL
1217 ($i, $c_i) = which_both($mask);
1219 This works just as "which", but the complement of $i will be in $c_i.
1221 pdl> $x = sequence(10); p $x
1223 [0 1 2 3 4 5 6 7 8 9]
1224 pdl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
1225 [5 6 7 8 9]
1226 [0 1 2 3 4]
1227 which_both processes bad values. It will set the bad-value flag of all
1229 output piddles if the flag is set for any of the input piddles.
1230 where
1232 Use a mask to select values from one or more data PDLs
1233 "where" accepts one or more data piddles and a mask piddle. It returns
1235 a list of output piddles, corresponding to the input data piddles.
1236 Each output piddle is a 1-dimensional list of values in its
1237 corresponding data piddle. The values are drawn from locations where
1238 the mask is nonzero.
1239 The output PDLs are still connected to the original data PDLs, for the
1241 purpose of dataflow.
1242 "where" combines the functionality of "which" and index into a single
1244 operation.
1245 BUGS:
1247 While "where" works OK for most N-dimensional cases, it does not thread
1249 properly over (for example) the (N+1)th dimension in data that is
1250 compared to an N-dimensional mask. Use "whereND" for that.
1251 $i = $x->where($x+5 > 0); # $i contains those elements of $x
1253 # where mask ($x+5 > 0) is 1
1254 $i .= -5; # Set those elements (of $x) to -5. Together, these
1255 # commands clamp $x to a maximum of -5.
1256 It is also possible to use the same mask for several piddles with the
1258 same call:
1259 ($i,$j,$k) = where($x,$y,$z, $x+5>0);
1261 Note: $i is always 1-D, even if $x is >1-D.
1263 WARNING: The first argument (the values) and the second argument (the
1265 mask) currently have to have the exact same dimensions (or horrible
1266 things happen). You *cannot* thread over a smaller mask, for example.
1267 whereND
1269 "where" with support for ND masks and threading
1270 "whereND" accepts one or more data piddles and a mask piddle. It
1272 returns a list of output piddles, corresponding to the input data
1273 piddles. The values are drawn from locations where the mask is
1274 nonzero.
1275 "whereND" differs from "where" in that the mask dimensionality is
1277 preserved which allows for proper threading of the selection operation
1278 over higher dimensions.
1279 As with "where" the output PDLs are still connected to the original
1281 data PDLs, for the purpose of dataflow.
1282 $sdata = whereND $data, $mask
1284 ($s1, $s2, ..., $sn) = whereND $d1, $d2, ..., $dn, $mask
1285 where
1287 $data is M dimensional
1289 $mask is N < M dimensional
1290 dims($data) 1..N == dims($mask) 1..N
1291 with threading over N+1 to M dimensions
1292 $data = sequence(4,3,2); # example data array
1294 $mask4 = (random(4)>0.5); # example 1-D mask array, has $n4 true values
1295 $mask43 = (random(4,3)>0.5); # example 2-D mask array, has $n43 true values
1296 $sdat4 = whereND $data, $mask4; # $sdat4 is a [$n4,3,2] pdl
1297 $sdat43 = whereND $data, $mask43; # $sdat43 is a [$n43,2] pdl
1298 Just as with "where", you can use the returned value in an assignment.
1300 That means that both of these examples are valid:
1301 # Used to create a new slice stored in $sdat4:
1303 $sdat4 = $data->whereND($mask4);
1304 $sdat4 .= 0;
1305 # Used in lvalue context:
1306 $data->whereND($mask4) .= 0;
1307 whichND
1309 Return the coordinates of non-zero values in a mask.
1310 WhichND returns the N-dimensional coordinates of each nonzero value in
1312 a mask PDL with any number of dimensions. The returned values arrive
1313 as an array-of-vectors suitable for use in indexND or range.
1314 $coords = whichND($mask);
1316 returns a PDL containing the coordinates of the elements that are non-
1318 zero in $mask, suitable for use in indexND. The 0th dimension contains
1319 the full coordinate listing of each point; the 1st dimension lists all
1320 the points. For example, if $mask has rank 4 and 100 matching
1321 elements, then $coords has dimension 4x100.
1322 If no such elements exist, then whichND returns a structured empty PDL:
1324 an Nx0 PDL that contains no values (but matches, threading-wise, with
1325 the vectors that would be produced if such elements existed).
1326 DEPRECATED BEHAVIOR IN LIST CONTEXT:
1328 whichND once delivered different values in list context than in scalar
1330 context, for historical reasons. In list context, it returned the
1331 coordinates transposed, as a collection of 1-PDLs (one per dimension)
1332 in a list. This usage is deprecated in PDL 2.4.10, and will cause a
1333 warning to be issued every time it is encountered. To avoid the
1334 warning, you can set the global variable "$PDL::whichND" to 's' to get
1335 scalar behavior in all contexts, or to 'l' to get list behavior in list
1336 context.
1337 In later versions of PDL, the deprecated behavior will disappear.
1339 Deprecated list context whichND expressions can be replaced with:
1340 @list = $x->whichND->mv(0,-1)->dog;
1342 SEE ALSO:
1344 "which" finds coordinates of nonzero values in a 1-D mask.
1346 "where" extracts values from a data PDL that are associated with
1348 nonzero values in a mask PDL.
1349 pdl> $s=sequence(10,10,3,4)
1351 pdl> ($x, $y, $z, $w)=whichND($s == 203); p $x, $y, $z, $w
1352 [3] [0] [2] [0]
1353 pdl> print $s->at(list(cat($x,$y,$z,$w)))
1354 203
1355 setops
1357 Implements simple set operations like union and intersection
1358 Usage: $set = setops($x, <OPERATOR>, $y);
1360 The operator can be "OR", "XOR" or "AND". This is then applied to $x
1362 viewed as a set and $y viewed as a set. Set theory says that a set may
1363 not have two or more identical elements, but setops takes care of this
1364 for you, so "$x=pdl(1,1,2)" is OK. The functioning is as follows:
1365 "OR"
1367 The resulting vector will contain the elements that are either in
1368 $x or in $y or both. This is the union in set operation terms
1369 "XOR"
1371 The resulting vector will contain the elements that are either in
1372 $x or $y, but not in both. This is
1373 Union($x, $y) - Intersection($x, $y)
1375 in set operation terms.
1377 "AND"
1379 The resulting vector will contain the intersection of $x and $y, so
1380 the elements that are in both $x and $y. Note that for convenience
1381 this operation is also aliased to "intersect".
1382 It should be emphasized that these routines are used when one or both
1384 of the sets $x, $y are hard to calculate or that you get from a
1385 separate subroutine.
1386 Finally IDL users might be familiar with Craig Markwardt's
1388 "cmset_op.pro" routine which has inspired this routine although it was
1389 written independently However the present routine has a few less
1390 options (but see the examples)
1391 You will very often use these functions on an index vector, so that is
1393 what we will show here. We will in fact something slightly silly. First
1394 we will find all squares that are also cubes below 10000.
1395 Create a sequence vector:
1397 pdl> $x = sequence(10000)
1399 Find all odd and even elements:
1401 pdl> ($even, $odd) = which_both( ($x % 2) == 0)
1403 Find all squares
1405 pdl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
1407 Find all cubes (being careful with roundoff error!)
1409 pdl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
1411 Then find all squares that are cubes:
1413 pdl> $both = setops($squares, 'AND', $cubes)
1415 And print these (assumes that "PDL::NiceSlice" is loaded!)
1417 pdl> p $x($both)
1419 [0 1 64 729 4096]
1420 Then find all numbers that are either cubes or squares, but not both:
1422 pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)
1424 pdl> p $cube_xor_square->nelem()
1426 112
1427 So there are a total of 112 of these!
1429 Finally find all odd squares:
1431 pdl> $odd_squares = setops($squares, 'AND', $odd)
1433 Another common occurrence is to want to get all objects that are in $x
1435 and in the complement of $y. But it is almost always best to create the
1436 complement explicitly since the universe that both are taken from is
1437 not known. Thus use "which_both" if possible to keep track of
1438 complements.
1439 If this is impossible the best approach is to make a temporary:
1441 This creates an index vector the size of the universe of the sets and
1443 set all elements in $y to 0
1444 pdl> $tmp = ones($n_universe); $tmp($y) .= 0;
1446 This then finds the complement of $y
1448 pdl> $C_b = which($tmp == 1);
1450 and this does the final selection:
1452 pdl> $set = setops($x, 'AND', $C_b)
1454 intersect
1456 Calculate the intersection of two piddles
1457 Usage: $set = intersect($x, $y);
1459 This routine is merely a simple interface to "setops". See that for
1461 more information
1462 Find all numbers less that 100 that are of the form 2*y and 3*x
1464 pdl> $x=sequence(100)
1466 pdl> $factor2 = which( ($x % 2) == 0)
1467 pdl> $factor3 = which( ($x % 3) == 0)
1468 pdl> $ii=intersect($factor2, $factor3)
1469 pdl> p $x($ii)
1470 [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
1471
1473 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
1474 Contributions by Christian Soeller (c.soeller@auckland.ac.nz), Karl
1475 Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
1476 (deforest@boulder.swri.edu) and Jarle Brinchmann (jarle@astro.up.pt)
1477 All rights reserved. There is no warranty. You are allowed to
1478 redistribute this software / documentation under certain conditions.
1479 For details, see the file COPYING in the PDL distribution. If this file
1480 is separated from the PDL distribution, the copyright notice should be
1481 included in the file.
1482 Updated for CPAN viewing compatibility by David Mertens.
1484perl v5.32.1 2021-02-15 Primitive(3)