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 "N**4".
174 The reason for having this routine is that you do not need to have the
176 same thread-dimensions for "tmp" as for the other arguments, which in
177 case of large numbers of matrices makes this much more memory-
178 efficient.
179 It is hoped that things like this could be taken care of as a kind of
181 closures at some point.
182 inner2t processes bad values. It will set the bad-value flag of all
184 output piddles if the flag is set for any of the input piddles.
185 crossp
187 Signature: (a(tri=3); b(tri); [o] c(tri))
188 Cross product of two 3D vectors
190 After
192 $c = crossp $x, $y
194 the inner product "$c*$x" and "$c*$y" will be zero, i.e. $c is
196 orthogonal to $x and $y
197 crossp does not process bad values. It will set the bad-value flag of
199 all output piddles if the flag is set for any of the input piddles.
200 norm
202 Signature: (vec(n); [o] norm(n))
203 Normalises a vector to unit Euclidean length
205 norm processes bad values. It will set the bad-value flag of all
207 output piddles if the flag is set for any of the input piddles.
208 indadd
210 Signature: (a(); indx ind(); [o] sum(m))
211 Threaded Index Add: Add "a" to the "ind" element of "sum", i.e:
213 sum(ind) += a
215 Simple Example:
217 $x = 2;
219 $ind = 3;
220 $sum = zeroes(10);
221 indadd($x,$ind, $sum);
222 print $sum
223 #Result: ( 2 added to element 3 of $sum)
224 # [0 0 0 2 0 0 0 0 0 0]
225 Threaded Example:
227 $x = pdl( 1,2,3);
229 $ind = pdl( 1,4,6);
230 $sum = zeroes(10);
231 indadd($x,$ind, $sum);
232 print $sum."\n";
233 #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
234 # [0 1 0 0 2 0 3 0 0 0]
235 The routine barfs if any of the indices are bad.
237 conv1d
239 Signature: (a(m); kern(p); [o]b(m); int reflect)
240 1D convolution along first dimension
242 The m-th element of the discrete convolution of an input piddle $a of
244 size $M, and a kernel piddle $kern of size $P, is calculated as
245 n = ($P-1)/2
247 ====
248 \
249 ($a conv1d $kern)[m] = > $a_ext[m - n] * $kern[n]
250 /
251 ====
252 n = -($P-1)/2
253 where $a_ext is either the periodic (or reflected) extension of $a so
255 it is equal to $a on " 0..$M-1 " and equal to the corresponding
256 periodic/reflected image of $a outside that range.
257 $con = conv1d sequence(10), pdl(-1,0,1);
259 $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
261 By default, periodic boundary conditions are assumed (i.e. wrap
263 around). Alternatively, you can request reflective boundary conditions
264 using the "Boundary" option:
265 {Boundary => 'reflect'} # case in 'reflect' doesn't matter
267 The convolution is performed along the first dimension. To apply it
269 across another dimension use the slicing routines, e.g.
270 $y = $x->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
272 This function is useful for threaded filtering of 1D signals.
274 Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv
276 WARNING: "conv1d" processes bad values in its inputs as the numeric
278 value of "$pdl->badvalue" so it is not recommended for processing pdls
279 with bad values in them unless special care is taken.
280 conv1d ignores the bad-value flag of the input piddles. It will set
282 the bad-value flag of all output piddles if the flag is set for any of
283 the input piddles.
284 in
286 Signature: (a(); b(n); [o] c())
287 test if a is in the set of values b
289 $goodmsk = $labels->in($goodlabels);
291 print pdl(3,1,4,6,2)->in(pdl(2,3,3));
292 [1 0 0 0 1]
293 "in" is akin to the is an element of of set theory. In principle, PDL
295 threading could be used to achieve its functionality by using a
296 construct like
297 $msk = ($labels->dummy(0) == $goodlabels)->orover;
299 However, "in" doesn't create a (potentially large) intermediate and is
301 generally faster.
302 in does not process bad values. It will set the bad-value flag of all
304 output piddles if the flag is set for any of the input piddles.
305 uniq
307 return all unique elements of a piddle
308 The unique elements are returned in ascending order.
310 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniq
312 [-1 0 2 4 6] # 0 is returned 2nd (sorted order)
313 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniq
315 [-1 2 4 6 nan] # NaN value is returned at end
316 Note: The returned pdl is 1D; any structure of the input piddle is
318 lost. "NaN" values are never compare equal to any other values, even
319 themselves. As a result, they are always unique. "uniq" returns the
320 NaN values at the end of the result piddle. This follows the Matlab
321 usage.
322 See uniqind if you need the indices of the unique elements rather than
324 the values.
325 Bad values are not considered unique by uniq and are ignored.
327 $x=sequence(10);
329 $x=$x->setbadif($x%3);
330 print $x->uniq;
331 [0 3 6 9]
332 uniqind
334 Return the indices of all unique elements of a piddle The order is in
335 the order of the values to be consistent with uniq. "NaN" values never
336 compare equal with any other value and so are always unique. This
337 follows the Matlab usage.
338 PDL> p pdl(2,2,2,4,0,-1,6,6)->uniqind
340 [5 4 1 3 6] # the 0 at index 4 is returned 2nd, but...
341 PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniqind
343 [5 1 3 6 4] # ...the NaN at index 4 is returned at end
344 Note: The returned pdl is 1D; any structure of the input piddle is
346 lost.
347 See uniq if you want the unique values instead of the indices.
349 Bad values are not considered unique by uniqind and are ignored.
351 uniqvec
353 Return all unique vectors out of a collection
354 NOTE: If any vectors in the input piddle have NaN values
356 they are returned at the end of the non-NaN ones. This is
357 because, by definition, NaN values never compare equal with
358 any other value.
359 NOTE: The current implementation does not sort the vectors
361 containing NaN values.
362 The unique vectors are returned in lexicographically sorted ascending
364 order. The 0th dimension of the input PDL is treated as a dimensional
365 index within each vector, and the 1st and any higher dimensions are
366 taken to run across vectors. The return value is always 2D; any
367 structure of the input PDL (beyond using the 0th dimension for vector
368 index) is lost.
369 See also uniq for a unique list of scalars; and qsortvec for sorting a
371 list of vectors lexicographcally.
372 If a vector contains all bad values, it is ignored as in uniq. If some
374 of the values are good, it is treated as a normal vector. For example,
375 [1 2 BAD] and [BAD 2 3] could be returned, but [BAD BAD BAD] could not.
376 Vectors containing BAD values will be returned after any non-NaN and
377 non-BAD containing vectors, followed by the NaN vectors.
378 hclip
380 Signature: (a(); b(); [o] c())
381 clip (threshold) $a by $b ($b is upper bound)
383 hclip processes bad values. It will set the bad-value flag of all
385 output piddles if the flag is set for any of the input piddles.
386 lclip
388 Signature: (a(); b(); [o] c())
389 clip (threshold) $a by $b ($b is lower bound)
391 lclip processes bad values. It will set the bad-value flag of all
393 output piddles if the flag is set for any of the input piddles.
394 clip
396 Clip (threshold) a piddle by (optional) upper or lower bounds.
397 $y = $x->clip(0,3);
399 $c = $x->clip(undef, $x);
400 clip handles bad values since it is just a wrapper around hclip and
402 lclip.
403 clip
405 Signature: (a(); l(); h(); [o] c())
406 info not available
408 clip processes bad values. It will set the bad-value flag of all
410 output piddles if the flag is set for any of the input piddles.
411 wtstat
413 Signature: (a(n); wt(n); avg(); [o]b(); int deg)
414 Weighted statistical moment of given degree
416 This calculates a weighted statistic over the vector "a". The formula
418 is
419 b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
421 Bad values are ignored in any calculation; $b will only have its bad
423 flag set if the output contains any bad data.
424 statsover
426 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())
427 Calculate useful statistics over a dimension of a piddle
429 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($piddle, $weights);
431 This utility function calculates various useful quantities of a piddle.
433 These are:
434 · the mean:
436 MEAN = sum (x)/ N
438 with "N" being the number of elements in x
440 · the population RMS deviation from the mean:
442 PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
444 The population deviation is the best-estimate of the deviation of
446 the population from which a sample is drawn.
447 · the median
449 The median is the 50th percentile data value. Median is found by
451 medover, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
452 · the minimum
454 · the maximum
456 · the average absolute deviation:
458 AADEV = sum( abs(x-mean(x)) )/N
460 · RMS deviation from the mean:
462 RMS = sqrt(sum( (x-mean(x))^2 )/N)
464 (also known as the root-mean-square deviation, or the square root of
466 the variance)
467 This operator is a projection operator so the calculation will take
469 place over the final dimension. Thus if the input is N-dimensional each
470 returned value will be N-1 dimensional, to calculate the statistics for
471 the entire piddle either use "clump(-1)" directly on the piddle or call
472 "stats".
473 Bad values are simply ignored in the calculation, effectively reducing
475 the sample size. If all data are bad then the output data are marked
476 bad.
477 stats
479 Calculates useful statistics on a piddle
480 ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($piddle,[$weights]);
482 This utility calculates all the most useful quantities in one call. It
484 works the same way as "statsover", except that the quantities are
485 calculated considering the entire input PDL as a single sample, rather
486 than as a collection of rows. See "statsover" for definitions of the
487 returned quantities.
488 Bad values are handled; if all input values are bad, then all of the
490 output values are flagged bad.
491 histogram
493 Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
494 Calculates a histogram for given stepsize and minimum.
496 $h = histogram($data, $step, $min, $numbins);
498 $hist = zeroes $numbins; # Put histogram in existing piddle.
499 histogram($data, $hist, $step, $min, $numbins);
500 The histogram will contain $numbins bins starting from $min, each $step
502 wide. The value in each bin is the number of values in $data that lie
503 within the bin limits.
504 Data below the lower limit is put in the first bin, and data above the
506 upper limit is put in the last bin.
507 The output is reset in a different threadloop so that you can take a
509 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
510 For a higher-level interface, see hist.
512 pdl> p histogram(pdl(1,1,2),1,0,3)
514 [0 2 1]
515 histogram processes bad values. It will set the bad-value flag of all
517 output piddles if the flag is set for any of the input piddles.
518 whistogram
520 Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
521 Calculates a histogram from weighted data for given stepsize and
523 minimum.
524 $h = whistogram($data, $weights, $step, $min, $numbins);
526 $hist = zeroes $numbins; # Put histogram in existing piddle.
527 whistogram($data, $weights, $hist, $step, $min, $numbins);
528 The histogram will contain $numbins bins starting from $min, each $step
530 wide. The value in each bin is the sum of the values in $weights that
531 correspond to values in $data that lie within the bin limits.
532 Data below the lower limit is put in the first bin, and data above the
534 upper limit is put in the last bin.
535 The output is reset in a different threadloop so that you can take a
537 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
538 pdl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
540 [0 0.2 0.5 0]
541 whistogram processes bad values. It will set the bad-value flag of all
543 output piddles if the flag is set for any of the input piddles.
544 histogram2d
546 Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
547 double stepb; double minb; int mbsize => mb;)
548 Calculates a 2d histogram.
550 $h = histogram2d($datax, $datay, $stepx, $minx,
552 $nbinx, $stepy, $miny, $nbiny);
553 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
554 histogram2d($datax, $datay, $hist, $stepx, $minx,
555 $nbinx, $stepy, $miny, $nbiny);
556 The histogram will contain $nbinx x $nbiny bins, with the lower limits
558 of the first one at "($minx, $miny)", and with bin size "($stepx,
559 $stepy)". The value in each bin is the number of values in $datax and
560 $datay that lie within the bin limits.
561 Data below the lower limit is put in the first bin, and data above the
563 upper limit is put in the last bin.
564 pdl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
566 [
567 [0 0 0]
568 [0 2 2]
569 [0 1 0]
570 ]
571 histogram2d processes bad values. It will set the bad-value flag of
573 all output piddles if the flag is set for any of the input piddles.
574 whistogram2d
576 Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
577 double stepb; double minb; int mbsize => mb;)
578 Calculates a 2d histogram from weighted data.
580 $h = whistogram2d($datax, $datay, $weights,
582 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
583 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
584 whistogram2d($datax, $datay, $weights, $hist,
585 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
586 The histogram will contain $nbinx x $nbiny bins, with the lower limits
588 of the first one at "($minx, $miny)", and with bin size "($stepx,
589 $stepy)". The value in each bin is the sum of the values in $weights
590 that correspond to values in $datax and $datay that lie within the bin
591 limits.
592 Data below the lower limit is put in the first bin, and data above the
594 upper limit is put in the last bin.
595 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)
597 [
598 [ 0 0 0]
599 [ 0 0.5 0.9]
600 [ 0 0.1 0]
601 ]
602 whistogram2d processes bad values. It will set the bad-value flag of
604 all output piddles if the flag is set for any of the input piddles.
605 fibonacci
607 Signature: ([o]x(n))
608 Constructor - a vector with Fibonacci's sequence
610 fibonacci does not process bad values. It will set the bad-value flag
612 of all output piddles if the flag is set for any of the input piddles.
613 append
615 Signature: (a(n); b(m); [o] c(mn))
616 append two piddles by concatenating along their first dimensions
618 $x = ones(2,4,7);
620 $y = sequence 5;
621 $c = $x->append($y); # size of $c is now (7,4,7) (a jumbo-piddle ;)
622 "append" appends two piddles along their first dimensions. The rest of
624 the dimensions must be compatible in the threading sense. The resulting
625 size of the first dimension is the sum of the sizes of the first
626 dimensions of the two argument piddles - i.e. "n + m".
627 Similar functions include glue (below), which can append more than two
629 piddles along an arbitrary dimension, and cat, which can append more
630 than two piddles that all have the same sized dimensions.
631 append does not process bad values. It will set the bad-value flag of
633 all output piddles if the flag is set for any of the input piddles.
634 glue
636 $c = $x->glue(<dim>,$y,...)
637 Glue two or more PDLs together along an arbitrary dimension (N-D
639 append).
640 Sticks $x, $y, and all following arguments together along the specified
642 dimension. All other dimensions must be compatible in the threading
643 sense.
644 Glue is permissive, in the sense that every PDL is treated as having an
646 infinite number of trivial dimensions of order 1 -- so "$x->glue(3,$y)"
647 works, even if $x and $y are only one dimensional.
648 If one of the PDLs has no elements, it is ignored. Likewise, if one of
650 them is actually the undefined value, it is treated as if it had no
651 elements.
652 If the first parameter is a defined perl scalar rather than a pdl, then
654 it is taken as a dimension along which to glue everything else, so you
655 can say "$cube = PDL::glue(3,@image_list);" if you like.
656 "glue" is implemented in pdl, using a combination of xchg and append.
658 It should probably be updated (one day) to a pure PP function.
659 Similar functions include append (above), which appends only two
661 piddles along their first dimension, and cat, which can append more
662 than two piddles that all have the same sized dimensions.
663 axisvalues
665 Signature: ([o,nc]a(n))
666 Internal routine
668 "axisvalues" is the internal primitive that implements axisvals and
670 alters its argument.
671 axisvalues does not process bad values. It will set the bad-value flag
673 of all output piddles if the flag is set for any of the input piddles.
674 random
676 Constructor which returns piddle of random numbers
677 $x = random([type], $nx, $ny, $nz,...);
679 $x = random $y;
680 etc (see zeroes).
682 This is the uniform distribution between 0 and 1 (assumedly excluding 1
684 itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
685 specify dimensions, types or give a template.
686 You can use the perl function srand to seed the random generator. For
688 further details consult Perl's srand documentation.
689 randsym
691 Constructor which returns piddle of random numbers
692 $x = randsym([type], $nx, $ny, $nz,...);
694 $x = randsym $y;
695 etc (see zeroes).
697 This is the uniform distribution between 0 and 1 (excluding both 0 and
699 1, cf random). The arguments are the same as "zeroes" (q.v.) - i.e. one
700 can specify dimensions, types or give a template.
701 You can use the perl function srand to seed the random generator. For
703 further details consult Perl's srand documentation.
704 grandom
706 Constructor which returns piddle of Gaussian random numbers
707 $x = grandom([type], $nx, $ny, $nz,...);
709 $x = grandom $y;
710 etc (see zeroes).
712 This is generated using the math library routine "ndtri".
714 Mean = 0, Stddev = 1
716 You can use the perl function srand to seed the random generator. For
718 further details consult Perl's srand documentation.
719 vsearch
721 Signature: ( vals(); xs(n); [o] indx(); [\%options] )
722 Efficiently search for values in a sorted piddle, returning indices.
724 $idx = vsearch( $vals, $x, [\%options] );
726 vsearch( $vals, $x, $idx, [\%options ] );
727 vsearch performs a binary search in the ordered piddle $x, for the
729 values from $vals piddle, returning indices into $x. What is a
730 "match", and the meaning of the returned indices, are determined by the
731 options.
732 The "mode" option indicates which method of searching to use, and may
734 be one of:
735 "sample"
737 invoke vsearch_sample, returning indices appropriate for sampling
738 within a distribution.
739 "insert_leftmost"
741 invoke vsearch_insert_leftmost, returning the left-most possible
742 insertion point which still leaves the piddle sorted.
743 "insert_rightmost"
745 invoke vsearch_insert_rightmost, returning the right-most possible
746 insertion point which still leaves the piddle sorted.
747 "match"
749 invoke vsearch_match, returning the index of a matching element,
750 else -(insertion point + 1)
751 "bin_inclusive"
753 invoke vsearch_bin_inclusive, returning an index appropriate for
754 binning on a grid where the left bin edges are inclusive of the
755 bin. See below for further explanation of the bin.
756 "bin_exclusive"
758 invoke vsearch_bin_exclusive, returning an index appropriate for
759 binning on a grid where the left bin edges are exclusive of the
760 bin. See below for further explanation of the bin.
761 The default value of "mode" is "sample".
763 use PDL;
765 my @modes = qw( sample insert_leftmost insert_rightmost match
767 bin_inclusive bin_exclusive );
768 # Generate a sequence of 3 zeros, 3 ones, ..., 3 fours.
770 my $x = zeroes(3,5)->yvals->flat;
771 for my $mode ( @modes ) {
773 # if the value is in $x
774 my $contained = 2;
775 my $idx_contained = vsearch( $contained, $x, { mode => $mode } );
776 my $x_contained = $x->copy;
777 $x_contained->slice( $idx_contained ) .= 9;
778 # if the value is not in $x
780 my $not_contained = 1.5;
781 my $idx_not_contained = vsearch( $not_contained, $x, { mode => $mode } );
782 my $x_not_contained = $x->copy;
783 $x_not_contained->slice( $idx_not_contained ) .= 9;
784 print sprintf("%-23s%30s\n", '$x', $x);
786 print sprintf("%-23s%30s\n", "$mode ($contained)", $x_contained);
787 print sprintf("%-23s%30s\n\n", "$mode ($not_contained)", $x_not_contained);
788 }
789 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
791 # sample (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
792 # sample (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
793 #
794 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
795 # insert_leftmost (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
796 # insert_leftmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
797 #
798 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
799 # insert_rightmost (2) [0 0 0 1 1 1 2 2 2 9 3 3 4 4 4]
800 # insert_rightmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
801 #
802 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
803 # match (2) [0 0 0 1 1 1 2 9 2 3 3 3 4 4 4]
804 # match (1.5) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
805 #
806 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
807 # bin_inclusive (2) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
808 # bin_inclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
809 #
810 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
811 # bin_exclusive (2) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
812 # bin_exclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
813 Also see vsearch_sample, vsearch_insert_leftmost,
815 vsearch_insert_rightmost, vsearch_match, vsearch_bin_inclusive, and
816 vsearch_bin_exclusive
817 vsearch_sample
819 Signature: (vals(); x(n); indx [o]idx())
820 Search for values in a sorted array, return index appropriate for
822 sampling from a distribution
823 $idx = vsearch_sample($vals, $x);
825 $x must be sorted, but may be in decreasing or increasing order.
827 vsearch_sample returns an index I for each value V of $vals appropriate
829 for sampling $vals
830 I has the following properties:
832 · if $x is sorted in increasing order
834 V <= x[0] : I = 0
836 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
837 x[-1] < V : I = $x->nelem -1
838 · if $x is sorted in decreasing order
840 V > x[0] : I = 0
842 x[0] >= V > x[-1] : I s.t. x[I] >= V > x[I+1]
843 x[-1] >= V : I = $x->nelem - 1
844 If all elements of $x are equal, I = $x->nelem - 1.
846 If $x contains duplicated elements, I is the index of the leftmost (by
848 position in array) duplicate if V matches.
849 This function is useful e.g. when you have a list of probabilities for
851 events and want to generate indices to events:
852 $x = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
854 $y = random 20;
855 $c = vsearch_sample($y, $x); # Now, $c will have the appropriate distr.
856 It is possible to use the cumusumover function to obtain cumulative
858 probabilities from absolute probabilities.
859 needs major (?) work to handles bad values
861 vsearch_insert_leftmost
863 Signature: (vals(); x(n); indx [o]idx())
864 Determine the insertion point for values in a sorted array, inserting
866 before duplicates.
867 $idx = vsearch_insert_leftmost($vals, $x);
869 $x must be sorted, but may be in decreasing or increasing order.
871 vsearch_insert_leftmost returns an index I for each value V of $vals
873 equal to the leftmost position (by index in array) within $x that V may
874 be inserted and still maintain the order in $x.
875 Insertion at index I involves shifting elements I and higher of $x to
877 the right by one and setting the now empty element at index I to V.
878 I has the following properties:
880 · if $x is sorted in increasing order
882 V <= x[0] : I = 0
884 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
885 x[-1] < V : I = $x->nelem
886 · if $x is sorted in decreasing order
888 V > x[0] : I = -1
890 x[0] >= V >= x[-1] : I s.t. x[I] >= V > x[I+1]
891 x[-1] >= V : I = $x->nelem -1
892 If all elements of $x are equal,
894 i = 0
896 If $x contains duplicated elements, I is the index of the leftmost (by
898 index in array) duplicate if V matches.
899 needs major (?) work to handles bad values
901 vsearch_insert_rightmost
903 Signature: (vals(); x(n); indx [o]idx())
904 Determine the insertion point for values in a sorted array, inserting
906 after duplicates.
907 $idx = vsearch_insert_rightmost($vals, $x);
909 $x must be sorted, but may be in decreasing or increasing order.
911 vsearch_insert_rightmost returns an index I for each value V of $vals
913 equal to the rightmost position (by index in array) within $x that V
914 may be inserted and still maintain the order in $x.
915 Insertion at index I involves shifting elements I and higher of $x to
917 the right by one and setting the now empty element at index I to V.
918 I has the following properties:
920 · if $x is sorted in increasing order
922 V < x[0] : I = 0
924 x[0] <= V < x[-1] : I s.t. x[I-1] <= V < x[I]
925 x[-1] <= V : I = $x->nelem
926 · if $x is sorted in decreasing order
928 V >= x[0] : I = -1
930 x[0] > V >= x[-1] : I s.t. x[I] >= V > x[I+1]
931 x[-1] > V : I = $x->nelem -1
932 If all elements of $x are equal,
934 i = $x->nelem - 1
936 If $x contains duplicated elements, I is the index of the leftmost (by
938 index in array) duplicate if V matches.
939 needs major (?) work to handles bad values
941 vsearch_match
943 Signature: (vals(); x(n); indx [o]idx())
944 Match values against a sorted array.
946 $idx = vsearch_match($vals, $x);
948 $x must be sorted, but may be in decreasing or increasing order.
950 vsearch_match returns an index I for each value V of $vals. If V
952 matches an element in $x, I is the index of that element, otherwise it
953 is -( insertion_point + 1 ), where insertion_point is an index in $x
954 where V may be inserted while maintaining the order in $x. If $x has
955 duplicated values, I may refer to any of them.
956 needs major (?) work to handles bad values
958 vsearch_bin_inclusive
960 Signature: (vals(); x(n); indx [o]idx())
961 Determine the index for values in a sorted array of bins, lower bound
963 inclusive.
964 $idx = vsearch_bin_inclusive($vals, $x);
966 $x must be sorted, but may be in decreasing or increasing order.
968 $x represents the edges of contiguous bins, with the first and last
970 elements representing the outer edges of the outer bins, and the inner
971 elements the shared bin edges.
972 The lower bound of a bin is inclusive to the bin, its outer bound is
974 exclusive to it. vsearch_bin_inclusive returns an index I for each
975 value V of $vals
976 I has the following properties:
978 · if $x is sorted in increasing order
980 V < x[0] : I = -1
982 x[0] <= V < x[-1] : I s.t. x[I] <= V < x[I+1]
983 x[-1] <= V : I = $x->nelem - 1
984 · if $x is sorted in decreasing order
986 V >= x[0] : I = 0
988 x[0] > V >= x[-1] : I s.t. x[I+1] > V >= x[I]
989 x[-1] > V : I = $x->nelem
990 If all elements of $x are equal,
992 i = $x->nelem - 1
994 If $x contains duplicated elements, I is the index of the righmost (by
996 index in array) duplicate if V matches.
997 needs major (?) work to handles bad values
999 vsearch_bin_exclusive
1001 Signature: (vals(); x(n); indx [o]idx())
1002 Determine the index for values in a sorted array of bins, lower bound
1004 exclusive.
1005 $idx = vsearch_bin_exclusive($vals, $x);
1007 $x must be sorted, but may be in decreasing or increasing order.
1009 $x represents the edges of contiguous bins, with the first and last
1011 elements representing the outer edges of the outer bins, and the inner
1012 elements the shared bin edges.
1013 The lower bound of a bin is exclusive to the bin, its upper bound is
1015 inclusive to it. vsearch_bin_exclusive returns an index I for each
1016 value V of $vals.
1017 I has the following properties:
1019 · if $x is sorted in increasing order
1021 V <= x[0] : I = -1
1023 x[0] < V <= x[-1] : I s.t. x[I] < V <= x[I+1]
1024 x[-1] < V : I = $x->nelem - 1
1025 · if $x is sorted in decreasing order
1027 V > x[0] : I = 0
1029 x[0] >= V > x[-1] : I s.t. x[I-1] >= V > x[I]
1030 x[-1] >= V : I = $x->nelem
1031 If all elements of $x are equal,
1033 i = $x->nelem - 1
1035 If $x contains duplicated elements, I is the index of the righmost (by
1037 index in array) duplicate if V matches.
1038 needs major (?) work to handles bad values
1040 interpolate
1042 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
1043 routine for 1D linear interpolation
1045 ( $yi, $err ) = interpolate($xi, $x, $y)
1047 Given a set of points "($x,$y)", use linear interpolation to find the
1049 values $yi at a set of points $xi.
1050 "interpolate" uses a binary search to find the suspects, er...,
1052 interpolation indices and therefore abscissas (ie $x) have to be
1053 strictly ordered (increasing or decreasing). For interpolation at lots
1054 of closely spaced abscissas an approach that uses the last index found
1055 as a start for the next search can be faster (compare Numerical Recipes
1056 "hunt" routine). Feel free to implement that on top of the binary
1057 search if you like. For out of bounds values it just does a linear
1058 extrapolation and sets the corresponding element of $err to 1, which is
1059 otherwise 0.
1060 See also interpol, which uses the same routine, differing only in the
1062 handling of extrapolation - an error message is printed rather than
1063 returning an error piddle.
1064 needs major (?) work to handles bad values
1066 interpol
1068 Signature: (xi(); x(n); y(n); [o] yi())
1069 routine for 1D linear interpolation
1071 $yi = interpol($xi, $x, $y)
1073 "interpol" uses the same search method as interpolate, hence $x must be
1075 strictly ordered (either increasing or decreasing). The difference
1076 occurs in the handling of out-of-bounds values; here an error message
1077 is printed.
1078 interpND
1080 Interpolate values from an N-D piddle, with switchable method
1081 $source = 10*xvals(10,10) + yvals(10,10);
1083 $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
1084 print $source->interpND( $index );
1085 InterpND acts like indexND, collapsing $index by lookup into $source;
1087 but it does interpolation rather than direct sampling. The
1088 interpolation method and boundary condition are switchable via an
1089 options hash.
1090 By default, linear or sample interpolation is used, with constant value
1092 outside the boundaries of the source pdl. No dataflow occurs, because
1093 in general the output is computed rather than indexed.
1094 All the interpolation methods treat the pixels as value-centered, so
1096 the "sample" method will return "$a->(0)" for coordinate values on the
1097 set [-0.5,0.5), and all methods will return "$a->(1)" for a coordinate
1098 value of exactly 1.
1099 Recognized options:
1101 method
1103 Values can be:
1104 · 0, s, sample, Sample (default for integer source types)
1106 The nearest value is taken. Pixels are regarded as centered on
1108 their respective integer coordinates (no offset from the linear
1109 case).
1110 · 1, l, linear, Linear (default for floating point source types)
1112 The values are N-linearly interpolated from an N-dimensional cube
1114 of size 2.
1115 · 3, c, cube, cubic, Cubic
1117 The values are interpolated using a local cubic fit to the data.
1119 The fit is constrained to match the original data and its
1120 derivative at the data points. The second derivative of the fit
1121 is not continuous at the data points. Multidimensional datasets
1122 are interpolated by the successive-collapse method.
1123 (Note that the constraint on the first derivative causes a small
1125 amount of ringing around sudden features such as step functions).
1126 · f, fft, fourier, Fourier
1128 The source is Fourier transformed, and the interpolated values
1130 are explicitly calculated from the coefficients. The boundary
1131 condition option is ignored -- periodic boundaries are imposed.
1132 If you pass in the option "fft", and it is a list (ARRAY) ref,
1134 then it is a stash for the magnitude and phase of the source FFT.
1135 If the list has two elements then they are taken as already
1136 computed; otherwise they are calculated and put in the stash.
1137 b, bound, boundary, Boundary
1139 This option is passed unmodified into indexND, which is used as the
1140 indexing engine for the interpolation. Some current allowed values
1141 are 'extend', 'periodic', 'truncate', and 'mirror' (default is
1142 'truncate').
1143 bad
1145 contains the fill value used for 'truncate' boundary. (default 0)
1146 fft
1148 An array ref whose associated list is used to stash the FFT of the
1149 source data, for the FFT method.
1150 one2nd
1152 Converts a one dimensional index piddle to a set of ND coordinates
1153 @coords=one2nd($x, $indices)
1155 returns an array of piddles containing the ND indexes corresponding to
1157 the one dimensional list indices. The indices are assumed to correspond
1158 to array $x clumped using "clump(-1)". This routine is used in the old
1159 vector form of whichND, but is useful on its own occasionally.
1160 Returned piddles have the indx datatype. $indices can have values
1162 larger than "$x->nelem" but negative values in $indices will not give
1163 the answer you expect.
1164 pdl> $x=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$x->clump(-1)
1166 pdl> $maxind=maximum_ind($c); p $maxind;
1167 6
1168 pdl> print one2nd($x, maximum_ind($c))
1169 0 1 1
1170 pdl> p $x->at(0,1,1)
1171 3
1172 which
1174 Signature: (mask(n); indx [o] inds(m))
1175 Returns indices of non-zero values from a 1-D PDL
1177 $i = which($mask);
1179 returns a pdl with indices for all those elements that are nonzero in
1181 the mask. Note that the returned indices will be 1D. If you feed in a
1182 multidimensional mask, it will be flattened before the indices are
1183 calculated. See also whichND for multidimensional masks.
1184 If you want to index into the original mask or a similar piddle with
1186 output from "which", remember to flatten it before calling index:
1187 $data = random 5, 5;
1189 $idx = which $data > 0.5; # $idx is now 1D
1190 $bigsum = $data->flat->index($idx)->sum; # flatten before indexing
1191 Compare also where for similar functionality.
1193 SEE ALSO:
1195 which_both returns separately the indices of both zero and nonzero
1197 values in the mask.
1198 where returns associated values from a data PDL, rather than indices
1200 into the mask PDL.
1201 whichND returns N-D indices into a multidimensional PDL.
1203 pdl> $x = sequence(10); p $x
1205 [0 1 2 3 4 5 6 7 8 9]
1206 pdl> $indx = which($x>6); p $indx
1207 [7 8 9]
1208 which processes bad values. It will set the bad-value flag of all
1210 output piddles if the flag is set for any of the input piddles.
1211 which_both
1213 Signature: (mask(n); indx [o] inds(m); indx [o]notinds(q))
1214 Returns indices of zero and nonzero values in a mask PDL
1216 ($i, $c_i) = which_both($mask);
1218 This works just as which, but the complement of $i will be in $c_i.
1220 pdl> $x = sequence(10); p $x
1222 [0 1 2 3 4 5 6 7 8 9]
1223 pdl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
1224 [5 6 7 8 9]
1225 [0 1 2 3 4]
1226 which_both processes bad values. It will set the bad-value flag of all
1228 output piddles if the flag is set for any of the input piddles.
1229 where
1231 Use a mask to select values from one or more data PDLs
1232 "where" accepts one or more data piddles and a mask piddle. It returns
1234 a list of output piddles, corresponding to the input data piddles.
1235 Each output piddle is a 1-dimensional list of values in its
1236 corresponding data piddle. The values are drawn from locations where
1237 the mask is nonzero.
1238 The output PDLs are still connected to the original data PDLs, for the
1240 purpose of dataflow.
1241 "where" combines the functionality of which and index into a single
1243 operation.
1244 BUGS:
1246 While "where" works OK for most N-dimensional cases, it does not thread
1248 properly over (for example) the (N+1)th dimension in data that is
1249 compared to an N-dimensional mask. Use "whereND" for that.
1250 $i = $x->where($x+5 > 0); # $i contains those elements of $x
1252 # where mask ($x+5 > 0) is 1
1253 $i .= -5; # Set those elements (of $x) to -5. Together, these
1254 # commands clamp $x to a maximum of -5.
1255 It is also possible to use the same mask for several piddles with the
1257 same call:
1258 ($i,$j,$k) = where($x,$y,$z, $x+5>0);
1260 Note: $i is always 1-D, even if $x is >1-D.
1262 WARNING: The first argument (the values) and the second argument (the
1264 mask) currently have to have the exact same dimensions (or horrible
1265 things happen). You *cannot* thread over a smaller mask, for example.
1266 whereND
1268 "where" with support for ND masks and threading
1269 "whereND" accepts one or more data piddles and a mask piddle. It
1271 returns a list of output piddles, corresponding to the input data
1272 piddles. The values are drawn from locations where the mask is
1273 nonzero.
1274 "whereND" differs from "where" in that the mask dimensionality is
1276 preserved which allows for proper threading of the selection operation
1277 over higher dimensions.
1278 As with "where" the output PDLs are still connected to the original
1280 data PDLs, for the purpose of dataflow.
1281 $sdata = whereND $data, $mask
1283 ($s1, $s2, ..., $sn) = whereND $d1, $d2, ..., $dn, $mask
1284 where
1286 $data is M dimensional
1288 $mask is N < M dimensional
1289 dims($data) 1..N == dims($mask) 1..N
1290 with threading over N+1 to M dimensions
1291 $data = sequence(4,3,2); # example data array
1293 $mask4 = (random(4)>0.5); # example 1-D mask array, has $n4 true values
1294 $mask43 = (random(4,3)>0.5); # example 2-D mask array, has $n43 true values
1295 $sdat4 = whereND $data, $mask4; # $sdat4 is a [$n4,3,2] pdl
1296 $sdat43 = whereND $data, $mask43; # $sdat43 is a [$n43,2] pdl
1297 Just as with "where", you can use the returned value in an assignment.
1299 That means that both of these examples are valid:
1300 # Used to create a new slice stored in $sdat4:
1302 $sdat4 = $data->whereND($mask4);
1303 $sdat4 .= 0;
1304 # Used in lvalue context:
1305 $data->whereND($mask4) .= 0;
1306 whichND
1308 Return the coordinates of non-zero values in a mask.
1309 WhichND returns the N-dimensional coordinates of each nonzero value in
1311 a mask PDL with any number of dimensions. The returned values arrive
1312 as an array-of-vectors suitable for use in indexND or range.
1313 $coords = whichND($mask);
1315 returns a PDL containing the coordinates of the elements that are non-
1317 zero in $mask, suitable for use in indexND. The 0th dimension contains
1318 the full coordinate listing of each point; the 1st dimension lists all
1319 the points. For example, if $mask has rank 4 and 100 matching
1320 elements, then $coords has dimension 4x100.
1321 If no such elements exist, then whichND returns a structured empty PDL:
1323 an Nx0 PDL that contains no values (but matches, threading-wise, with
1324 the vectors that would be produced if such elements existed).
1325 DEPRECATED BEHAVIOR IN LIST CONTEXT:
1327 whichND once delivered different values in list context than in scalar
1329 context, for historical reasons. In list context, it returned the
1330 coordinates transposed, as a collection of 1-PDLs (one per dimension)
1331 in a list. This usage is deprecated in PDL 2.4.10, and will cause a
1332 warning to be issued every time it is encountered. To avoid the
1333 warning, you can set the global variable "$PDL::whichND" to 's' to get
1334 scalar behavior in all contexts, or to 'l' to get list behavior in list
1335 context.
1336 In later versions of PDL, the deprecated behavior will disappear.
1338 Deprecated list context whichND expressions can be replaced with:
1339 @list = $x->whichND->mv(0,-1)->dog;
1341 SEE ALSO:
1343 which finds coordinates of nonzero values in a 1-D mask.
1345 where extracts values from a data PDL that are associated with nonzero
1347 values in a mask PDL.
1348 pdl> $s=sequence(10,10,3,4)
1350 pdl> ($x, $y, $z, $w)=whichND($s == 203); p $x, $y, $z, $w
1351 [3] [0] [2] [0]
1352 pdl> print $s->at(list(cat($x,$y,$z,$w)))
1353 203
1354 setops
1356 Implements simple set operations like union and intersection
1357 Usage: $set = setops($x, <OPERATOR>, $y);
1359 The operator can be "OR", "XOR" or "AND". This is then applied to $x
1361 viewed as a set and $y viewed as a set. Set theory says that a set may
1362 not have two or more identical elements, but setops takes care of this
1363 for you, so "$x=pdl(1,1,2)" is OK. The functioning is as follows:
1364 "OR"
1366 The resulting vector will contain the elements that are either in
1367 $x or in $y or both. This is the union in set operation terms
1368 "XOR"
1370 The resulting vector will contain the elements that are either in
1371 $x or $y, but not in both. This is
1372 Union($x, $y) - Intersection($x, $y)
1374 in set operation terms.
1376 "AND"
1378 The resulting vector will contain the intersection of $x and $y, so
1379 the elements that are in both $x and $y. Note that for convenience
1380 this operation is also aliased to intersect.
1381 It should be emphasized that these routines are used when one or both
1383 of the sets $x, $y are hard to calculate or that you get from a
1384 separate subroutine.
1385 Finally IDL users might be familiar with Craig Markwardt's
1387 "cmset_op.pro" routine which has inspired this routine although it was
1388 written independently However the present routine has a few less
1389 options (but see the examples)
1390 You will very often use these functions on an index vector, so that is
1392 what we will show here. We will in fact something slightly silly. First
1393 we will find all squares that are also cubes below 10000.
1394 Create a sequence vector:
1396 pdl> $x = sequence(10000)
1398 Find all odd and even elements:
1400 pdl> ($even, $odd) = which_both( ($x % 2) == 0)
1402 Find all squares
1404 pdl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
1406 Find all cubes (being careful with roundoff error!)
1408 pdl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
1410 Then find all squares that are cubes:
1412 pdl> $both = setops($squares, 'AND', $cubes)
1414 And print these (assumes that "PDL::NiceSlice" is loaded!)
1416 pdl> p $x($both)
1418 [0 1 64 729 4096]
1419 Then find all numbers that are either cubes or squares, but not both:
1421 pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)
1423 pdl> p $cube_xor_square->nelem()
1425 112
1426 So there are a total of 112 of these!
1428 Finally find all odd squares:
1430 pdl> $odd_squares = setops($squares, 'AND', $odd)
1432 Another common occurrence is to want to get all objects that are in $x
1434 and in the complement of $y. But it is almost always best to create the
1435 complement explicitly since the universe that both are taken from is
1436 not known. Thus use which_both if possible to keep track of
1437 complements.
1438 If this is impossible the best approach is to make a temporary:
1440 This creates an index vector the size of the universe of the sets and
1442 set all elements in $y to 0
1443 pdl> $tmp = ones($n_universe); $tmp($y) .= 0;
1445 This then finds the complement of $y
1447 pdl> $C_b = which($tmp == 1);
1449 and this does the final selection:
1451 pdl> $set = setops($x, 'AND', $C_b)
1453 intersect
1455 Calculate the intersection of two piddles
1456 Usage: $set = intersect($x, $y);
1458 This routine is merely a simple interface to setops. See that for more
1460 information
1461 Find all numbers less that 100 that are of the form 2*y and 3*x
1463 pdl> $x=sequence(100)
1465 pdl> $factor2 = which( ($x % 2) == 0)
1466 pdl> $factor3 = which( ($x % 3) == 0)
1467 pdl> $ii=intersect($factor2, $factor3)
1468 pdl> p $x($ii)
1469 [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
1470
1472 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
1473 Contributions by Christian Soeller (c.soeller@auckland.ac.nz), Karl
1474 Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
1475 (deforest@boulder.swri.edu) and Jarle Brinchmann (jarle@astro.up.pt)
1476 All rights reserved. There is no warranty. You are allowed to
1477 redistribute this software / documentation under certain conditions.
1478 For details, see the file COPYING in the PDL distribution. If this file
1479 is separated from the PDL distribution, the copyright notice should be
1480 included in the file.
1481 Updated for CPAN viewing compatibility by David Mertens.
1483perl v5.32.0 2020-09-17 Primitive(3)