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 ndarrays if the flag is set for any of the input ndarrays.
45 x
47 Signature: (a(i,z), b(x,i),[o]c(x,z))
48 Matrix multiplication
50 PDL overloads the "x" operator (normally the repeat operator) for
52 matrix multiplication. The number of columns (size of the 0 dimension)
53 in the left-hand argument must normally equal the number of rows (size
54 of the 1 dimension) in the right-hand argument.
55 Row vectors are represented as (N x 1) two-dimensional PDLs, or you may
57 be sloppy and use a one-dimensional PDL. Column vectors are
58 represented as (1 x N) two-dimensional PDLs.
59 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 ndarrays. It will set
129 the bad-value flag of all output ndarrays if the flag is set for any of
130 the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarray $a of
245 size $M, and a kernel ndarray $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 ndarrays. It will set
283 the bad-value flag of all output ndarrays if the flag is set for any of
284 the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
306 uniq
308 return all unique elements of an ndarray
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 ndarray 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 ndarray. 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 an ndarray 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 ndarray 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 ndarray 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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarrays if the flag is set for any of the input ndarrays.
395 clip
397 Clip (threshold) an ndarray 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 ndarrays if the flag is set for any of the input ndarrays.
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 an ndarray
430 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($ndarray, $weights);
432 This utility function calculates various useful quantities of an
434 ndarray. 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 ndarray either use "clump(-1)" directly on the ndarray or
473 call "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 an ndarray
481 ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($ndarray,[$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 ndarray.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarray.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarray.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarray.
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 ndarrays if the flag is set for any of the input ndarrays.
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 ndarrays if the flag is set for any of the input
614 ndarrays.
615 append
617 Signature: (a(n); b(m); [o] c(mn))
618 append two ndarrays by concatenating along their first dimensions
620 $x = ones(2,4,7);
622 $y = sequence 5;
623 $c = $x->append($y); # size of $c is now (7,4,7) (a jumbo-ndarray ;)
624 "append" appends two ndarrays along their first dimensions. The rest of
626 the dimensions must be compatible in the threading sense. The resulting
627 size of the first dimension is the sum of the sizes of the first
628 dimensions of the two argument ndarrays - i.e. "n + m".
629 Similar functions include "glue" (below), which can append more than
631 two ndarrays along an arbitrary dimension, and cat, which can append
632 more than two ndarrays that all have the same sized dimensions.
633 append does not process bad values. It will set the bad-value flag of
635 all output ndarrays if the flag is set for any of the input ndarrays.
636 glue
638 $c = $x->glue(<dim>,$y,...)
639 Glue two or more PDLs together along an arbitrary dimension (N-D
641 "append").
642 Sticks $x, $y, and all following arguments together along the specified
644 dimension. All other dimensions must be compatible in the threading
645 sense.
646 Glue is permissive, in the sense that every PDL is treated as having an
648 infinite number of trivial dimensions of order 1 -- so "$x->glue(3,$y)"
649 works, even if $x and $y are only one dimensional.
650 If one of the PDLs has no elements, it is ignored. Likewise, if one of
652 them is actually the undefined value, it is treated as if it had no
653 elements.
654 If the first parameter is a defined perl scalar rather than a pdl, then
656 it is taken as a dimension along which to glue everything else, so you
657 can say "$cube = PDL::glue(3,@image_list);" if you like.
658 "glue" is implemented in pdl, using a combination of xchg and "append".
660 It should probably be updated (one day) to a pure PP function.
661 Similar functions include "append" (above), which appends only two
663 ndarrays along their first dimension, and cat, which can append more
664 than two ndarrays that all have the same sized dimensions.
665 axisvalues
667 Signature: ([o,nc]a(n))
668 Internal routine
670 "axisvalues" is the internal primitive that implements axisvals and
672 alters its argument.
673 axisvalues does not process bad values. It will set the bad-value flag
675 of all output ndarrays if the flag is set for any of the input
676 ndarrays.
677 random
679 Constructor which returns ndarray of random numbers
680 $x = random([type], $nx, $ny, $nz,...);
682 $x = random $y;
683 etc (see zeroes).
685 This is the uniform distribution between 0 and 1 (assumedly excluding 1
687 itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
688 specify dimensions, types or give a template.
689 You can use the perl function srand to seed the random generator. For
691 further details consult Perl's srand documentation.
692 randsym
694 Constructor which returns ndarray of random numbers
695 $x = randsym([type], $nx, $ny, $nz,...);
697 $x = randsym $y;
698 etc (see zeroes).
700 This is the uniform distribution between 0 and 1 (excluding both 0 and
702 1, cf "random"). The arguments are the same as "zeroes" (q.v.) - i.e.
703 one can specify dimensions, types or give a template.
704 You can use the perl function srand to seed the random generator. For
706 further details consult Perl's srand documentation.
707 grandom
709 Constructor which returns ndarray of Gaussian random numbers
710 $x = grandom([type], $nx, $ny, $nz,...);
712 $x = grandom $y;
713 etc (see zeroes).
715 This is generated using the math library routine "ndtri".
717 Mean = 0, Stddev = 1
719 You can use the perl function srand to seed the random generator. For
721 further details consult Perl's srand documentation.
722 vsearch
724 Signature: ( vals(); xs(n); [o] indx(); [\%options] )
725 Efficiently search for values in a sorted ndarray, returning indices.
727 $idx = vsearch( $vals, $x, [\%options] );
729 vsearch( $vals, $x, $idx, [\%options ] );
730 vsearch performs a binary search in the ordered ndarray $x, for the
732 values from $vals ndarray, returning indices into $x. What is a
733 "match", and the meaning of the returned indices, are determined by the
734 options.
735 The "mode" option indicates which method of searching to use, and may
737 be one of:
738 "sample"
740 invoke vsearch_sample, returning indices appropriate for sampling
741 within a distribution.
742 "insert_leftmost"
744 invoke vsearch_insert_leftmost, returning the left-most possible
745 insertion point which still leaves the ndarray sorted.
746 "insert_rightmost"
748 invoke vsearch_insert_rightmost, returning the right-most possible
749 insertion point which still leaves the ndarray sorted.
750 "match"
752 invoke vsearch_match, returning the index of a matching element,
753 else -(insertion point + 1)
754 "bin_inclusive"
756 invoke vsearch_bin_inclusive, returning an index appropriate for
757 binning on a grid where the left bin edges are inclusive of the
758 bin. See below for further explanation of the bin.
759 "bin_exclusive"
761 invoke vsearch_bin_exclusive, returning an index appropriate for
762 binning on a grid where the left bin edges are exclusive of the
763 bin. See below for further explanation of the bin.
764 The default value of "mode" is "sample".
766 use PDL;
768 my @modes = qw( sample insert_leftmost insert_rightmost match
770 bin_inclusive bin_exclusive );
771 # Generate a sequence of 3 zeros, 3 ones, ..., 3 fours.
773 my $x = zeroes(3,5)->yvals->flat;
774 for my $mode ( @modes ) {
776 # if the value is in $x
777 my $contained = 2;
778 my $idx_contained = vsearch( $contained, $x, { mode => $mode } );
779 my $x_contained = $x->copy;
780 $x_contained->slice( $idx_contained ) .= 9;
781 # if the value is not in $x
783 my $not_contained = 1.5;
784 my $idx_not_contained = vsearch( $not_contained, $x, { mode => $mode } );
785 my $x_not_contained = $x->copy;
786 $x_not_contained->slice( $idx_not_contained ) .= 9;
787 print sprintf("%-23s%30s\n", '$x', $x);
789 print sprintf("%-23s%30s\n", "$mode ($contained)", $x_contained);
790 print sprintf("%-23s%30s\n\n", "$mode ($not_contained)", $x_not_contained);
791 }
792 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
794 # sample (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
795 # sample (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
796 #
797 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
798 # insert_leftmost (2) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
799 # insert_leftmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
800 #
801 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
802 # insert_rightmost (2) [0 0 0 1 1 1 2 2 2 9 3 3 4 4 4]
803 # insert_rightmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
804 #
805 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
806 # match (2) [0 0 0 1 1 1 2 9 2 3 3 3 4 4 4]
807 # match (1.5) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
808 #
809 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
810 # bin_inclusive (2) [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
811 # bin_inclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
812 #
813 # $x [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
814 # bin_exclusive (2) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
815 # bin_exclusive (1.5) [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
816 Also see vsearch_sample, vsearch_insert_leftmost,
818 vsearch_insert_rightmost, vsearch_match, vsearch_bin_inclusive, and
819 vsearch_bin_exclusive
820 vsearch_sample
822 Signature: (vals(); x(n); indx [o]idx())
823 Search for values in a sorted array, return index appropriate for
825 sampling from a distribution
826 $idx = vsearch_sample($vals, $x);
828 $x must be sorted, but may be in decreasing or increasing order.
830 vsearch_sample returns an index I for each value V of $vals appropriate
832 for sampling $vals
833 I has the following properties:
835 • if $x is sorted in increasing order
837 V <= x[0] : I = 0
839 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
840 x[-1] < V : I = $x->nelem -1
841 • if $x is sorted in decreasing order
843 V > x[0] : I = 0
845 x[0] >= V > x[-1] : I s.t. x[I] >= V > x[I+1]
846 x[-1] >= V : I = $x->nelem - 1
847 If all elements of $x are equal, I = $x->nelem - 1.
849 If $x contains duplicated elements, I is the index of the leftmost (by
851 position in array) duplicate if V matches.
852 This function is useful e.g. when you have a list of probabilities for
854 events and want to generate indices to events:
855 $x = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
857 $y = random 20;
858 $c = vsearch_sample($y, $x); # Now, $c will have the appropriate distr.
859 It is possible to use the cumusumover function to obtain cumulative
861 probabilities from absolute probabilities.
862 needs major (?) work to handles bad values
864 vsearch_insert_leftmost
866 Signature: (vals(); x(n); indx [o]idx())
867 Determine the insertion point for values in a sorted array, inserting
869 before duplicates.
870 $idx = vsearch_insert_leftmost($vals, $x);
872 $x must be sorted, but may be in decreasing or increasing order.
874 vsearch_insert_leftmost returns an index I for each value V of $vals
876 equal to the leftmost position (by index in array) within $x that V may
877 be inserted and still maintain the order in $x.
878 Insertion at index I involves shifting elements I and higher of $x to
880 the right by one and setting the now empty element at index I to V.
881 I has the following properties:
883 • if $x is sorted in increasing order
885 V <= x[0] : I = 0
887 x[0] < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
888 x[-1] < V : I = $x->nelem
889 • if $x is sorted in decreasing order
891 V > x[0] : I = -1
893 x[0] >= V >= x[-1] : I s.t. x[I] >= V > x[I+1]
894 x[-1] >= V : I = $x->nelem -1
895 If all elements of $x are equal,
897 i = 0
899 If $x contains duplicated elements, I is the index of the leftmost (by
901 index in array) duplicate if V matches.
902 needs major (?) work to handles bad values
904 vsearch_insert_rightmost
906 Signature: (vals(); x(n); indx [o]idx())
907 Determine the insertion point for values in a sorted array, inserting
909 after duplicates.
910 $idx = vsearch_insert_rightmost($vals, $x);
912 $x must be sorted, but may be in decreasing or increasing order.
914 vsearch_insert_rightmost returns an index I for each value V of $vals
916 equal to the rightmost position (by index in array) within $x that V
917 may be inserted and still maintain the order in $x.
918 Insertion at index I involves shifting elements I and higher of $x to
920 the right by one and setting the now empty element at index I to V.
921 I has the following properties:
923 • if $x is sorted in increasing order
925 V < x[0] : I = 0
927 x[0] <= V < x[-1] : I s.t. x[I-1] <= V < x[I]
928 x[-1] <= V : I = $x->nelem
929 • if $x is sorted in decreasing order
931 V >= x[0] : I = -1
933 x[0] > V >= x[-1] : I s.t. x[I] >= V > x[I+1]
934 x[-1] > V : I = $x->nelem -1
935 If all elements of $x are equal,
937 i = $x->nelem - 1
939 If $x contains duplicated elements, I is the index of the leftmost (by
941 index in array) duplicate if V matches.
942 needs major (?) work to handles bad values
944 vsearch_match
946 Signature: (vals(); x(n); indx [o]idx())
947 Match values against a sorted array.
949 $idx = vsearch_match($vals, $x);
951 $x must be sorted, but may be in decreasing or increasing order.
953 vsearch_match returns an index I for each value V of $vals. If V
955 matches an element in $x, I is the index of that element, otherwise it
956 is -( insertion_point + 1 ), where insertion_point is an index in $x
957 where V may be inserted while maintaining the order in $x. If $x has
958 duplicated values, I may refer to any of them.
959 needs major (?) work to handles bad values
961 vsearch_bin_inclusive
963 Signature: (vals(); x(n); indx [o]idx())
964 Determine the index for values in a sorted array of bins, lower bound
966 inclusive.
967 $idx = vsearch_bin_inclusive($vals, $x);
969 $x must be sorted, but may be in decreasing or increasing order.
971 $x represents the edges of contiguous bins, with the first and last
973 elements representing the outer edges of the outer bins, and the inner
974 elements the shared bin edges.
975 The lower bound of a bin is inclusive to the bin, its outer bound is
977 exclusive to it. vsearch_bin_inclusive returns an index I for each
978 value V of $vals
979 I has the following properties:
981 • if $x is sorted in increasing order
983 V < x[0] : I = -1
985 x[0] <= V < x[-1] : I s.t. x[I] <= V < x[I+1]
986 x[-1] <= V : I = $x->nelem - 1
987 • if $x is sorted in decreasing order
989 V >= x[0] : I = 0
991 x[0] > V >= x[-1] : I s.t. x[I+1] > V >= x[I]
992 x[-1] > V : I = $x->nelem
993 If all elements of $x are equal,
995 i = $x->nelem - 1
997 If $x contains duplicated elements, I is the index of the righmost (by
999 index in array) duplicate if V matches.
1000 needs major (?) work to handles bad values
1002 vsearch_bin_exclusive
1004 Signature: (vals(); x(n); indx [o]idx())
1005 Determine the index for values in a sorted array of bins, lower bound
1007 exclusive.
1008 $idx = vsearch_bin_exclusive($vals, $x);
1010 $x must be sorted, but may be in decreasing or increasing order.
1012 $x represents the edges of contiguous bins, with the first and last
1014 elements representing the outer edges of the outer bins, and the inner
1015 elements the shared bin edges.
1016 The lower bound of a bin is exclusive to the bin, its upper bound is
1018 inclusive to it. vsearch_bin_exclusive returns an index I for each
1019 value V of $vals.
1020 I has the following properties:
1022 • if $x is sorted in increasing order
1024 V <= x[0] : I = -1
1026 x[0] < V <= x[-1] : I s.t. x[I] < V <= x[I+1]
1027 x[-1] < V : I = $x->nelem - 1
1028 • if $x is sorted in decreasing order
1030 V > x[0] : I = 0
1032 x[0] >= V > x[-1] : I s.t. x[I-1] >= V > x[I]
1033 x[-1] >= V : I = $x->nelem
1034 If all elements of $x are equal,
1036 i = $x->nelem - 1
1038 If $x contains duplicated elements, I is the index of the righmost (by
1040 index in array) duplicate if V matches.
1041 needs major (?) work to handles bad values
1043 interpolate
1045 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
1046 routine for 1D linear interpolation
1048 ( $yi, $err ) = interpolate($xi, $x, $y)
1050 Given a set of points "($x,$y)", use linear interpolation to find the
1052 values $yi at a set of points $xi.
1053 "interpolate" uses a binary search to find the suspects, er...,
1055 interpolation indices and therefore abscissas (ie $x) have to be
1056 strictly ordered (increasing or decreasing). For interpolation at lots
1057 of closely spaced abscissas an approach that uses the last index found
1058 as a start for the next search can be faster (compare Numerical Recipes
1059 "hunt" routine). Feel free to implement that on top of the binary
1060 search if you like. For out of bounds values it just does a linear
1061 extrapolation and sets the corresponding element of $err to 1, which is
1062 otherwise 0.
1063 See also "interpol", which uses the same routine, differing only in the
1065 handling of extrapolation - an error message is printed rather than
1066 returning an error ndarray.
1067 needs major (?) work to handles bad values
1069 interpol
1071 Signature: (xi(); x(n); y(n); [o] yi())
1072 routine for 1D linear interpolation
1074 $yi = interpol($xi, $x, $y)
1076 "interpol" uses the same search method as "interpolate", hence $x must
1078 be strictly ordered (either increasing or decreasing). The difference
1079 occurs in the handling of out-of-bounds values; here an error message
1080 is printed.
1081 interpND
1083 Interpolate values from an N-D ndarray, with switchable method
1084 $source = 10*xvals(10,10) + yvals(10,10);
1086 $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
1087 print $source->interpND( $index );
1088 InterpND acts like indexND, collapsing $index by lookup into $source;
1090 but it does interpolation rather than direct sampling. The
1091 interpolation method and boundary condition are switchable via an
1092 options hash.
1093 By default, linear or sample interpolation is used, with constant value
1095 outside the boundaries of the source pdl. No dataflow occurs, because
1096 in general the output is computed rather than indexed.
1097 All the interpolation methods treat the pixels as value-centered, so
1099 the "sample" method will return "$a->(0)" for coordinate values on the
1100 set [-0.5,0.5), and all methods will return "$a->(1)" for a coordinate
1101 value of exactly 1.
1102 Recognized options:
1104 method
1106 Values can be:
1107 • 0, s, sample, Sample (default for integer source types)
1109 The nearest value is taken. Pixels are regarded as centered on
1111 their respective integer coordinates (no offset from the linear
1112 case).
1113 • 1, l, linear, Linear (default for floating point source types)
1115 The values are N-linearly interpolated from an N-dimensional cube
1117 of size 2.
1118 • 3, c, cube, cubic, Cubic
1120 The values are interpolated using a local cubic fit to the data.
1122 The fit is constrained to match the original data and its
1123 derivative at the data points. The second derivative of the fit
1124 is not continuous at the data points. Multidimensional datasets
1125 are interpolated by the successive-collapse method.
1126 (Note that the constraint on the first derivative causes a small
1128 amount of ringing around sudden features such as step functions).
1129 • f, fft, fourier, Fourier
1131 The source is Fourier transformed, and the interpolated values
1133 are explicitly calculated from the coefficients. The boundary
1134 condition option is ignored -- periodic boundaries are imposed.
1135 If you pass in the option "fft", and it is a list (ARRAY) ref,
1137 then it is a stash for the magnitude and phase of the source FFT.
1138 If the list has two elements then they are taken as already
1139 computed; otherwise they are calculated and put in the stash.
1140 b, bound, boundary, Boundary
1142 This option is passed unmodified into indexND, which is used as the
1143 indexing engine for the interpolation. Some current allowed values
1144 are 'extend', 'periodic', 'truncate', and 'mirror' (default is
1145 'truncate').
1146 bad
1148 contains the fill value used for 'truncate' boundary. (default 0)
1149 fft
1151 An array ref whose associated list is used to stash the FFT of the
1152 source data, for the FFT method.
1153 one2nd
1155 Converts a one dimensional index ndarray to a set of ND coordinates
1156 @coords=one2nd($x, $indices)
1158 returns an array of ndarrays containing the ND indexes corresponding to
1160 the one dimensional list indices. The indices are assumed to correspond
1161 to array $x clumped using "clump(-1)". This routine is used in the old
1162 vector form of "whichND", but is useful on its own occasionally.
1163 Returned ndarrays have the indx datatype. $indices can have values
1165 larger than "$x->nelem" but negative values in $indices will not give
1166 the answer you expect.
1167 pdl> $x=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$x->clump(-1)
1169 pdl> $maxind=maximum_ind($c); p $maxind;
1170 6
1171 pdl> print one2nd($x, maximum_ind($c))
1172 0 1 1
1173 pdl> p $x->at(0,1,1)
1174 3
1175 which
1177 Signature: (mask(n); indx [o] inds(m))
1178 Returns indices of non-zero values from a 1-D PDL
1180 $i = which($mask);
1182 returns a pdl with indices for all those elements that are nonzero in
1184 the mask. Note that the returned indices will be 1D. If you feed in a
1185 multidimensional mask, it will be flattened before the indices are
1186 calculated. See also "whichND" for multidimensional masks.
1187 If you want to index into the original mask or a similar ndarray with
1189 output from "which", remember to flatten it before calling index:
1190 $data = random 5, 5;
1192 $idx = which $data > 0.5; # $idx is now 1D
1193 $bigsum = $data->flat->index($idx)->sum; # flatten before indexing
1194 Compare also "where" for similar functionality.
1196 SEE ALSO:
1198 "which_both" returns separately the indices of both zero and nonzero
1200 values in the mask.
1201 "where" returns associated values from a data PDL, rather than indices
1203 into the mask PDL.
1204 "whichND" returns N-D indices into a multidimensional PDL.
1206 pdl> $x = sequence(10); p $x
1208 [0 1 2 3 4 5 6 7 8 9]
1209 pdl> $indx = which($x>6); p $indx
1210 [7 8 9]
1211 which processes bad values. It will set the bad-value flag of all
1213 output ndarrays if the flag is set for any of the input ndarrays.
1214 which_both
1216 Signature: (mask(n); indx [o] inds(m); indx [o]notinds(q))
1217 Returns indices of zero and nonzero values in a mask PDL
1219 ($i, $c_i) = which_both($mask);
1221 This works just as "which", but the complement of $i will be in $c_i.
1223 pdl> $x = sequence(10); p $x
1225 [0 1 2 3 4 5 6 7 8 9]
1226 pdl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
1227 [5 6 7 8 9]
1228 [0 1 2 3 4]
1229 which_both processes bad values. It will set the bad-value flag of all
1231 output ndarrays if the flag is set for any of the input ndarrays.
1232 where
1234 Use a mask to select values from one or more data PDLs
1235 "where" accepts one or more data ndarrays and a mask ndarray. It
1237 returns a list of output ndarrays, corresponding to the input data
1238 ndarrays. Each output ndarray is a 1-dimensional list of values in its
1239 corresponding data ndarray. The values are drawn from locations where
1240 the mask is nonzero.
1241 The output PDLs are still connected to the original data PDLs, for the
1243 purpose of dataflow.
1244 "where" combines the functionality of "which" and index into a single
1246 operation.
1247 BUGS:
1249 While "where" works OK for most N-dimensional cases, it does not thread
1251 properly over (for example) the (N+1)th dimension in data that is
1252 compared to an N-dimensional mask. Use "whereND" for that.
1253 $i = $x->where($x+5 > 0); # $i contains those elements of $x
1255 # where mask ($x+5 > 0) is 1
1256 $i .= -5; # Set those elements (of $x) to -5. Together, these
1257 # commands clamp $x to a maximum of -5.
1258 It is also possible to use the same mask for several ndarrays with the
1260 same call:
1261 ($i,$j,$k) = where($x,$y,$z, $x+5>0);
1263 Note: $i is always 1-D, even if $x is >1-D.
1265 WARNING: The first argument (the values) and the second argument (the
1267 mask) currently have to have the exact same dimensions (or horrible
1268 things happen). You *cannot* thread over a smaller mask, for example.
1269 whereND
1271 "where" with support for ND masks and threading
1272 "whereND" accepts one or more data ndarrays and a mask ndarray. It
1274 returns a list of output ndarrays, corresponding to the input data
1275 ndarrays. The values are drawn from locations where the mask is
1276 nonzero.
1277 "whereND" differs from "where" in that the mask dimensionality is
1279 preserved which allows for proper threading of the selection operation
1280 over higher dimensions.
1281 As with "where" the output PDLs are still connected to the original
1283 data PDLs, for the purpose of dataflow.
1284 $sdata = whereND $data, $mask
1286 ($s1, $s2, ..., $sn) = whereND $d1, $d2, ..., $dn, $mask
1287 where
1289 $data is M dimensional
1291 $mask is N < M dimensional
1292 dims($data) 1..N == dims($mask) 1..N
1293 with threading over N+1 to M dimensions
1294 $data = sequence(4,3,2); # example data array
1296 $mask4 = (random(4)>0.5); # example 1-D mask array, has $n4 true values
1297 $mask43 = (random(4,3)>0.5); # example 2-D mask array, has $n43 true values
1298 $sdat4 = whereND $data, $mask4; # $sdat4 is a [$n4,3,2] pdl
1299 $sdat43 = whereND $data, $mask43; # $sdat43 is a [$n43,2] pdl
1300 Just as with "where", you can use the returned value in an assignment.
1302 That means that both of these examples are valid:
1303 # Used to create a new slice stored in $sdat4:
1305 $sdat4 = $data->whereND($mask4);
1306 $sdat4 .= 0;
1307 # Used in lvalue context:
1308 $data->whereND($mask4) .= 0;
1309 whichND
1311 Return the coordinates of non-zero values in a mask.
1312 WhichND returns the N-dimensional coordinates of each nonzero value in
1314 a mask PDL with any number of dimensions. The returned values arrive
1315 as an array-of-vectors suitable for use in indexND or range.
1316 $coords = whichND($mask);
1318 returns a PDL containing the coordinates of the elements that are non-
1320 zero in $mask, suitable for use in indexND. The 0th dimension contains
1321 the full coordinate listing of each point; the 1st dimension lists all
1322 the points. For example, if $mask has rank 4 and 100 matching
1323 elements, then $coords has dimension 4x100.
1324 If no such elements exist, then whichND returns a structured empty PDL:
1326 an Nx0 PDL that contains no values (but matches, threading-wise, with
1327 the vectors that would be produced if such elements existed).
1328 DEPRECATED BEHAVIOR IN LIST CONTEXT:
1330 whichND once delivered different values in list context than in scalar
1332 context, for historical reasons. In list context, it returned the
1333 coordinates transposed, as a collection of 1-PDLs (one per dimension)
1334 in a list. This usage is deprecated in PDL 2.4.10, and will cause a
1335 warning to be issued every time it is encountered. To avoid the
1336 warning, you can set the global variable "$PDL::whichND" to 's' to get
1337 scalar behavior in all contexts, or to 'l' to get list behavior in list
1338 context.
1339 In later versions of PDL, the deprecated behavior will disappear.
1341 Deprecated list context whichND expressions can be replaced with:
1342 @list = $x->whichND->mv(0,-1)->dog;
1344 SEE ALSO:
1346 "which" finds coordinates of nonzero values in a 1-D mask.
1348 "where" extracts values from a data PDL that are associated with
1350 nonzero values in a mask PDL.
1351 pdl> $s=sequence(10,10,3,4)
1353 pdl> ($x, $y, $z, $w)=whichND($s == 203); p $x, $y, $z, $w
1354 [3] [0] [2] [0]
1355 pdl> print $s->at(list(cat($x,$y,$z,$w)))
1356 203
1357 setops
1359 Implements simple set operations like union and intersection
1360 Usage: $set = setops($x, <OPERATOR>, $y);
1362 The operator can be "OR", "XOR" or "AND". This is then applied to $x
1364 viewed as a set and $y viewed as a set. Set theory says that a set may
1365 not have two or more identical elements, but setops takes care of this
1366 for you, so "$x=pdl(1,1,2)" is OK. The functioning is as follows:
1367 "OR"
1369 The resulting vector will contain the elements that are either in
1370 $x or in $y or both. This is the union in set operation terms
1371 "XOR"
1373 The resulting vector will contain the elements that are either in
1374 $x or $y, but not in both. This is
1375 Union($x, $y) - Intersection($x, $y)
1377 in set operation terms.
1379 "AND"
1381 The resulting vector will contain the intersection of $x and $y, so
1382 the elements that are in both $x and $y. Note that for convenience
1383 this operation is also aliased to "intersect".
1384 It should be emphasized that these routines are used when one or both
1386 of the sets $x, $y are hard to calculate or that you get from a
1387 separate subroutine.
1388 Finally IDL users might be familiar with Craig Markwardt's
1390 "cmset_op.pro" routine which has inspired this routine although it was
1391 written independently However the present routine has a few less
1392 options (but see the examples)
1393 You will very often use these functions on an index vector, so that is
1395 what we will show here. We will in fact something slightly silly. First
1396 we will find all squares that are also cubes below 10000.
1397 Create a sequence vector:
1399 pdl> $x = sequence(10000)
1401 Find all odd and even elements:
1403 pdl> ($even, $odd) = which_both( ($x % 2) == 0)
1405 Find all squares
1407 pdl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
1409 Find all cubes (being careful with roundoff error!)
1411 pdl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
1413 Then find all squares that are cubes:
1415 pdl> $both = setops($squares, 'AND', $cubes)
1417 And print these (assumes that "PDL::NiceSlice" is loaded!)
1419 pdl> p $x($both)
1421 [0 1 64 729 4096]
1422 Then find all numbers that are either cubes or squares, but not both:
1424 pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)
1426 pdl> p $cube_xor_square->nelem()
1428 112
1429 So there are a total of 112 of these!
1431 Finally find all odd squares:
1433 pdl> $odd_squares = setops($squares, 'AND', $odd)
1435 Another common occurrence is to want to get all objects that are in $x
1437 and in the complement of $y. But it is almost always best to create the
1438 complement explicitly since the universe that both are taken from is
1439 not known. Thus use "which_both" if possible to keep track of
1440 complements.
1441 If this is impossible the best approach is to make a temporary:
1443 This creates an index vector the size of the universe of the sets and
1445 set all elements in $y to 0
1446 pdl> $tmp = ones($n_universe); $tmp($y) .= 0;
1448 This then finds the complement of $y
1450 pdl> $C_b = which($tmp == 1);
1452 and this does the final selection:
1454 pdl> $set = setops($x, 'AND', $C_b)
1456 intersect
1458 Calculate the intersection of two ndarrays
1459 Usage: $set = intersect($x, $y);
1461 This routine is merely a simple interface to "setops". See that for
1463 more information
1464 Find all numbers less that 100 that are of the form 2*y and 3*x
1466 pdl> $x=sequence(100)
1468 pdl> $factor2 = which( ($x % 2) == 0)
1469 pdl> $factor3 = which( ($x % 3) == 0)
1470 pdl> $ii=intersect($factor2, $factor3)
1471 pdl> p $x($ii)
1472 [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
1473
1475 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
1476 Contributions by Christian Soeller (c.soeller@auckland.ac.nz), Karl
1477 Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
1478 (deforest@boulder.swri.edu) and Jarle Brinchmann (jarle@astro.up.pt)
1479 All rights reserved. There is no warranty. You are allowed to
1480 redistribute this software / documentation under certain conditions.
1481 For details, see the file COPYING in the PDL distribution. If this file
1482 is separated from the PDL distribution, the copyright notice should be
1483 included in the file.
1484 Updated for CPAN viewing compatibility by David Mertens.
1486perl v5.34.0 2021-08-16 Primitive(3)