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 use PDL::Primitive;
17
19 inner
20 Signature: (a(n); b(n); [o]c())
22 Inner product over one dimension
24 c = sum_i a_i * b_i
26 outer
28 Signature: (a(n); b(m); [o]c(n,m))
30 outer product over one dimension
32 Naturally, it is possible to achieve the effects of outer product sim‐
34 ply by threading over the ""*"" operator but this function is provided
35 for convenience.
36 x
38 Signature: (a(i,x), b(z,i),[o]c(x,z))
40 Matrix multiplication
42 PDL overloads the "x" operator (normally the repeat operator) for
44 matrix multiplication. The number of columns (size of the 0 dimension)
45 in the left-hand argument must normally equal the number of rows (size
46 of the 1 dimension) in the right-hand argument.
47 Row vectors are represented as (N x 1) two-dimensional PDLs, or you may
49 be sloppy and use a one-dimensional PDL. Column vectors are repre‐
50 sented as (1 x N) two-dimensional PDLs.
51 Threading occurs in the usual way, but as both the 0 and 1 dimension
53 (if present) are included in the operation, you must be sure that you
54 don't try to thread over either of those dims.
55 EXAMPLES
57 Here are some simple ways to define vectors and matrices:
59 perldl> $r = pdl(1,2); # A row vector
61 perldl> $c = pdl([[3],[4]]); # A column vector
62 perldl> $c = pdl(3,4)->(*1); # A column vector, using NiceSlice
63 perldl> $m = pdl([[1,2],[3,4]]); # A 2x2 matrix
64 Now that we have a few objects prepared, here is how to matrix-multiply
66 them:
67 perldl> print $r x $m # row x matrix = row
69 [
70 [ 7 10]
71 ]
72 perldl> print $m x $r # matrix x row = ERROR
74 PDL: Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1
75 perldl> print $m x $c # matrix x column = column
77 [
78 [ 5]
79 [11]
80 ]
81 perldl> print $m x 2 # Trivial case: scalar mult.
83 [
84 [2 4]
85 [6 8]
86 ]
87 perldl> print $r x $c # row x column = scalar
89 [
90 [11]
91 ]
92 perldl> print $c x $r # column x row = matrix
94 [
95 [3 6]
96 [4 8]
97 ]
98 INTERNALS
100 The mechanics of the multiplication are carried out by the matmult
102 method.
103 matmult
105 Signature: (a(i,x),b(z,i),[o]c(x,z))
107 Matrix multiplication
109 We peruse the inner product to define matrix multiplication via a
111 threaded inner product.
112 For usage, see x, a description of the overloaded 'x' operator
114 innerwt
116 Signature: (a(n); b(n); c(n); [o]d())
118 Weighted (i.e. triple) inner product
120 d = sum_i a(i) b(i) c(i)
122 inner2
124 Signature: (a(n); b(n,m); c(m); [o]d())
126 Inner product of two vectors and a matrix
128 d = sum_ij a(i) b(i,j) c(j)
130 Note that you should probably not thread over "a" and "c" since that
132 would be very wasteful. Instead, you should use a temporary for "b*c".
133 inner2d
135 Signature: (a(n,m); b(n,m); [o]c())
137 Inner product over 2 dimensions.
139 Equivalent to
141 $c = inner($a->clump(2), $b->clump(2))
143 inner2t
145 Signature: (a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k)))
147 Efficient Triple matrix product "a*b*c"
149 Efficiency comes from by using the temporary "tmp". This operation only
151 scales as "N**3" whereas threading using inner2 would scale as "N**4".
152 The reason for having this routine is that you do not need to have the
154 same thread-dimensions for "tmp" as for the other arguments, which in
155 case of large numbers of matrices makes this much more memory-effi‐
156 cient.
157 It is hoped that things like this could be taken care of as a kind of
159 closures at some point.
160 crossp
162 Signature: (a(tri=3); b(tri); [o] c(tri))
164 Cross product of two 3D vectors
166 After
168 $c = crossp $a, $b
170 the inner product "$c*$a" and "$c*$b" will be zero, i.e. $c is orthogo‐
172 nal to $a and $b
173 norm
175 Signature: (vec(n); [o] norm(n))
177 Normalises a vector to unit Euclidean length
179 indadd
181 Signature: (a(); int ind(); [o] sum(m))
183 Threaded Index Add: Add "a" to the "ind" element of "sum", i.e:
185 sum(ind) += a
187 Simple Example:
189 $a = 2;
191 $ind = 3;
192 $sum = zeroes(10);
193 indadd($a,$ind, $sum);
194 print $sum
195 #Result: ( 2 added to element 3 of $sum)
196 # [0 0 0 2 0 0 0 0 0 0]
197 Threaded Example:
199 $a = pdl( 1,2,3);
201 $ind = pdl( 1,4,6);
202 $sum = zeroes(10);
203 indadd($a,$ind, $sum);
204 print $sum."\n";
205 #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
206 # [0 1 0 0 2 0 3 0 0 0]
207 conv1d
209 Signature: (a(m); kern(p); [o]b(m); int reflect)
211 1d convolution along first dimension
213 $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
215 By default, periodic boundary conditions are assumed (i.e. wrap
217 around). Alternatively, you can request reflective boundary conditions
218 using the "Boundary" option:
219 {Boundary => 'reflect'} # case in 'reflect' doesn't matter
221 The convolution is performed along the first dimension. To apply it
223 across another dimension use the slicing routines, e.g.
224 $b = $a->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
226 This function is useful for threaded filtering of 1D signals.
228 Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv
230 in
232 Signature: (a(); b(n); [o] c())
234 test if a is in the set of values b
236 $goodmsk = $labels->in($goodlabels);
238 print pdl(4,3,1)->in(pdl(2,3,3));
239 [0 1 0]
240 "in" is akin to the is an element of of set theory. In priciple, PDL
242 threading could be used to achieve its functionality by using a con‐
243 struct like
244 $msk = ($labels->dummy(0) == $goodlabels)->orover;
246 However, "in" doesn't create a (potentially large) intermediate and is
248 generally faster.
249 uniq
251 return all unique elements of a piddle
253 The unique elements are returned in ascending order.
255 print pdl(2,2,2,4,0,-1,6,6)->uniq;
257 [-1 0 2 4 6]
258 Note: The returned pdl is 1D; any structure of the input piddle is
260 lost.
261 See uniqind if you need the indices of the unique elements rather than
263 the values.
264 uniqind
266 return the indices of all unique elements of a piddle The order is in
268 the order of the values to be consistent with uniq
269 print pdl(2,2,2,4,0,-1,6,6)->uniqind;
271 [5, 4, 1, 3, 6]
272 Note: The returned pdl is 1D; any structure of the input piddle is
274 lost.
275 See uniq if you want the unique values instead of the indices.
277 uniqvec
279 return all unique vectors out of a collection
281 The unique vectors are returned in lexicographically sorted ascending
283 order. The 0th dimension of the input PDL is treated as a dimensional
284 index within each vector, and the 1st and any higher dimensions are
285 taken to run across vectors. The return value is always 2D; any struc‐
286 ture of the input PDL (beyond using the 0th dimension for vector index)
287 is lost.
288 See also uniq for a uniqe list of scalars; and qsortvec for sorting a
290 list of vectors lexicographcally.
291 hclip
293 Signature: (a(); b(); [o] c())
295 clip (threshold) $a by $b ($b is upper bound)
297 lclip
299 Signature: (a(); b(); [o] c())
301 clip (threshold) $a by $b ($b is lower bound)
303 clip
305 Clip (threshold) a piddle by (optional) upper or lower bounds.
307 $b = $a->clip(0,3);
309 $c = $a->clip(undef, $x);
310 wtstat
312 Signature: (a(n); wt(n); avg(); [o]b(); int deg)
314 Weighted statistical moment of given degree
316 This calculates a weighted statistic over the vector "a". The formula
318 is
319 b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
321 statsover
323 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())
325 Calculate useful statistics over a dimension of a piddle
327 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($piddle, $weights);
329 This utility function calculates various useful quantities of a piddle.
331 These are:
332 * the mean:
334 MEAN = sum (x)/ N
335 with "N" being the number of elements in x
337 * RMS deviation from the mean:
339 RMS = sqrt(sum( (x-mean(x))^2 )/N)
340 (also known as the root-mean-square deviation, or the square root of
342 the variance)
343 * the median
345 The median is the 50th percentile data value. Median is found by
346 medover, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
347 * the minimum
349 * the maximum
350 * the absolute deviation:
351 ADEV = sqrt(sum( abs(x-mean(x)) )/N)
352 (This is also called the standard deviation)
354 * the population RMS deviation from the mean:
356 PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
357 The population deviation is the best-estimate of the deviation of
359 the population from which a sample is drawn.
360 This operator is a projection operator so the calculation will take
362 place over the final dimension. Thus if the input is N-dimensional each
363 returned value will be N-1 dimensional, to calculate the statistics for
364 the entire piddle either use "clump(-1)" directly on the piddle or call
365 "stats".
366 stats
368 Calculates useful statistics on a piddle
370 ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($piddle,[$weights]);
372 This utility calculates all the most useful quantities in one call. It
374 works the same way as "statsover", except that the quantities are cal‐
375 culated considering the entire input PDL as a single sample, rather
376 than as a collection of rows. See "statsover" for definitions of the
377 returned quantities.
378 histogram
380 Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
382 Calculates a histogram for given stepsize and minimum.
384 $h = histogram($data, $step, $min, $numbins);
386 $hist = zeroes $numbins; # Put histogram in existing piddle.
387 histogram($data, $hist, $step, $min, $numbins);
388 The histogram will contain $numbins bins starting from $min, each $step
390 wide. The value in each bin is the number of values in $data that lie
391 within the bin limits.
392 Data below the lower limit is put in the first bin, and data above the
394 upper limit is put in the last bin.
395 The output is reset in a different threadloop so that you can take a
397 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
398 Use hist instead for a high-level interface.
400 perldl> p histogram(pdl(1,1,2),1,0,3)
402 [0 2 1]
403 whistogram
405 Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
407 Calculates a histogram from weighted data for given stepsize and mini‐
409 mum.
410 $h = whistogram($data, $weights, $step, $min, $numbins);
412 $hist = zeroes $numbins; # Put histogram in existing piddle.
413 whistogram($data, $weights, $hist, $step, $min, $numbins);
414 The histogram will contain $numbins bins starting from $min, each $step
416 wide. The value in each bin is the sum of the values in $weights that
417 correspond to values in $data that lie within the bin limits.
418 Data below the lower limit is put in the first bin, and data above the
420 upper limit is put in the last bin.
421 The output is reset in a different threadloop so that you can take a
423 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
424 perldl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
426 [0 0.2 0.5 0]
427 histogram2d
429 Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
431 double stepb; double minb; int mbsize => mb;)
432 Calculates a 2d histogram.
434 $h = histogram2d($datax, $datay,
436 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
437 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
438 histogram2d($datax, $datay, $hist,
439 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
440 The histogram will contain $nbinx x $nbiny bins, with the lower limits
442 of the first one at "($minx, $miny)", and with bin size "($stepx,
443 $stepy)". The value in each bin is the number of values in $datax and
444 $datay that lie within the bin limits.
445 Data below the lower limit is put in the first bin, and data above the
447 upper limit is put in the last bin.
448 perldl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
450 [
451 [0 0 0]
452 [0 2 2]
453 [0 1 0]
454 ]
455 whistogram2d
457 Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
459 double stepb; double minb; int mbsize => mb;)
460 Calculates a 2d histogram from weighted data.
462 $h = whistogram2d($datax, $datay, $weights,
464 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
465 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
466 whistogram2d($datax, $datay, $weights, $hist,
467 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
468 The histogram will contain $nbinx x $nbiny bins, with the lower limits
470 of the first one at "($minx, $miny)", and with bin size "($stepx,
471 $stepy)". The value in each bin is the sum of the values in $weights
472 that correspond to values in $datax and $datay that lie within the bin
473 limits.
474 Data below the lower limit is put in the first bin, and data above the
476 upper limit is put in the last bin.
477 perldl> 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)
479 [
480 [ 0 0 0]
481 [ 0 0.5 0.9]
482 [ 0 0.1 0]
483 ]
484 fibonacci
486 Signature: ([o]x(n))
488 Constructor - a vector with Fibonacci's sequence
490 append
492 Signature: (a(n); b(m); [o] c(mn))
494 append two or more piddles by concatenating along their first dimen‐
496 sions
497 $a = ones(2,4,7);
499 $b = sequence 5;
500 $c = $a->append($b); # size of $c is now (7,4,7) (a jumbo-piddle ;)
501 "append" appends two piddles along their first dims. Rest of the dimen‐
503 sions must be compatible in the threading sense. Resulting size of
504 first dim is the sum of the sizes of the first dims of the two argument
505 piddles - ie "n + m".
506 glue
508 $c = $a->glue(<dim>,$b,...)
510 Glue two or more PDLs together along an arbitrary dimension (N-D
512 append).
513 Sticks $a, $b, and all following arguments together along the specified
515 dimension. All other dimensions must be compatible in the threading
516 sense.
517 Glue is permissive, in the sense that every PDL is treated as having an
519 infinite number of trivial dimensions of order 1 -- so "$a-"glue(3,$b)>
520 works, even if $a and $b are only one dimensional.
521 If one of the PDLs has no elements, it is ignored. Likewise, if one of
523 them is actually the undefined value, it is treated as if it had no
524 elements.
525 If the first parameter is a defined perl scalar rather than a pdl, then
527 it is taken as a dimension along which to glue everything else, so you
528 can say "$cube = PDL::glue(3,@image_list);" if you like.
529 "glue" is implemented in pdl, using a combination of xchg and append.
531 It should probably be updated (one day) to a pure PP function.
532 axisvalues
534 Signature: ([o,nc]a(n))
536 Internal routine
538 "axisvalues" is the internal primitive that implements axisvals and
540 alters its argument.
541 random
543 Constructor which returns piddle of random numbers
545 $a = random([type], $nx, $ny, $nz,...);
547 $a = random $b;
548 etc (see zeroes).
550 This is the uniform distribution between 0 and 1 (assumedly excluding 1
552 itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
553 specify dimensions, types or give a template.
554 You can use the perl function srand to seed the random generator. For
556 further details consult Perl's srand documentation.
557 randsym
559 Constructor which returns piddle of random numbers
561 $a = randsym([type], $nx, $ny, $nz,...);
563 $a = randsym $b;
564 etc (see zeroes).
566 This is the uniform distribution between 0 and 1 (excluding both 0 and
568 1, cf random). The arguments are the same as "zeroes" (q.v.) - i.e. one
569 can specify dimensions, types or give a template.
570 You can use the perl function srand to seed the random generator. For
572 further details consult Perl's srand documentation.
573 grandom
575 Constructor which returns piddle of Gaussian random numbers
577 $a = grandom([type], $nx, $ny, $nz,...);
579 $a = grandom $b;
580 etc (see zeroes).
582 This is generated using the math library routine "ndtri".
584 Mean = 0, Stddev = 1
586 You can use the perl function srand to seed the random generator. For
588 further details consult Perl's srand documentation.
589 vsearch
591 Signature: (i(); x(n); int [o]ip())
593 routine for searching 1D values i.e. step-function interpolation.
595 $inds = vsearch($vals, $xs);
597 Returns for each value of $vals the index of the least larger member of
599 $xs (which need to be in increasing order). If the value is larger than
600 any member of $xs, the index to the last element of $xs is returned.
601 This function is useful e.g. when you have a list of probabilities for
603 events and want to generate indices to events:
604 $a = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
606 $b = random 20;
607 $c = vsearch($b, $a); # Now, $c will have the appropriate distr.
608 It is possible to use the cumusumover function to obtain cumulative
610 probabilities from absolute probabilities.
611 interpolate
613 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
615 routine for 1D linear interpolation
617 ( $yi, $err ) = interpolate($xi, $x, $y)
619 Given a set of points "($x,$y)", use linear interpolation to find the
621 values $yi at a set of points $xi.
622 "interpolate" uses a binary search to find the suspects, er..., inter‐
624 polation indices and therefore abscissas (ie $x) have to be strictly
625 ordered (increasing or decreasing). For interpolation at lots of
626 closely spaced abscissas an approach that uses the last index found as
627 a start for the next search can be faster (compare Numerical Recipes
628 "hunt" routine). Feel free to implement that on top of the binary
629 search if you like. For out of bounds values it just does a linear
630 extrapolation and sets the corresponding element of $err to 1, which is
631 otherwise 0.
632 See also interpol, which uses the same routine, differing only in the
634 handling of extrapolation - an error message is printed rather than
635 returning an error piddle.
636 interpol
638 Signature: (xi(); x(n); y(n); [o] yi())
640 routine for 1D linear interpolation
642 $yi = interpol($xi, $x, $y)
644 "interpol" uses the same search method as interpolate, hence $x must be
646 strictly ordered (either increasing or decreasing). The difference
647 occurs in the handling of out-of-bounds values; here an error message
648 is printed.
649 interpND
651 Interpolate values from an N-D piddle, with switchable method
653 $source = 10*xvals(10,10) + yvals(10,10);
655 $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
656 print $source->interpND( $index );
657 InterpND acts like indexND, collapsing $index by lookup into $source;
659 but it does interpolation rather than direct sampling. The interpola‐
660 tion method and boundary condition are switchable via an options hash.
661 By default, linear or sample interpolation is used, with constant value
663 outside the boundaries of the source pdl. No dataflow occurs, because
664 in general the output is computed rather than indexed.
665 All the interpolation methods treat the pixels as value-centered, so
667 the "sample" method will return $a->(0) for coordinate values on the
668 set [-0.5,0.5), and all methods will return $a->(1) for a coordinate
669 value of exactly 1.
670 Recognized options:
672 method
674 Values can be:
675 * 0, s, sample, Sample (default for integer source types)
677 The nearest value is taken. Pixels are regarded as centered on
678 their respective integer coordinates (no offset from the linear
679 case).
680 * 1, l, linear, Linear (default for floating point source types)
682 The values are N-linearly interpolated from an N-dimensional cube
683 of size 2.
684 * 3, c, cube, cubic, Cubic
686 The values are interpolated using a local cubic fit to the data.
687 The fit is constrained to match the original data and its deriva‐
688 tive at the data points. The second derivative of the fit is not
689 continuous at the data points. Multidimensional datasets are
690 interpolated by the successive-collapse method.
691 (Note that the constraint on the first derivative causes a small
693 amount of ringing around sudden features such as step functions).
694 * f, fft, fourier, Fourier
696 The source is Fourier transformed, and the interpolated values
697 are explicitly calculated from the coefficients. The boundary
698 condition option is ignored -- periodic boundaries are imposed.
699 If you pass in the option "fft", and it is a list (ARRAY) ref,
701 then it is a stash for the magnitude and phase of the source FFT.
702 If the list has two elements then they are taken as already com‐
703 puted; otherwise they are calculated and put in the stash.
704 b, bound, boundary, Boundary
706 This option is passed unmodified into indexND, which is used as the
707 indexing engine for the interpolation. Some current allowed values
708 are 'extend', 'periodic', 'truncate', and 'mirror' (default is
709 'truncate').
710 bad
712 contains the fill value used for 'truncate' boundary. (default 0)
713 fft
715 An array ref whose associated list is used to stash the FFT of the
716 source data, for the FFT method.
717 one2nd
719 Converts a one dimensional index piddle to a set of ND coordinates
721 @coords=one2nd($a, $indices)
723 returns an array of piddles containing the ND indexes corresponding to
725 the one dimensional list indices. The indices are assumed to correspond
726 to array $a clumped using "clump(-1)". This routine is used in whichND,
727 but is useful on its own occasionally.
728 perldl> $a=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$a->clump(-1)
730 perldl> $maxind=maximum_ind($c); p $maxind;
731 6
732 perldl> print one2nd($a, maximum_ind($c))
733 0 1 1
734 perldl> p $a->at(0,1,1)
735 3
736 which
738 Signature: (mask(n); int [o] inds(m))
740 Returns indices of non-zero values from a 1-D PDL
742 $i = which($mask);
744 returns a pdl with indices for all those elements that are nonzero in
746 the mask. Note that the returned indices will be 1D. If you feed in a
747 multidimensional mask, it will be flattened before the indices are cal‐
748 culated. See also whichND for multidimensional masks.
749 If you want to index into the original mask or a similar piddle with
751 output from "which", remember to flatten it before calling index:
752 $data = random 5, 5;
754 $idx = which $data > 0.5; # $idx is now 1D
755 $bigsum = $data->flat->index($idx)->sum; # flatten before indexing
756 Compare also where for similar functionality.
758 SEE ALSO:
760 which_both returns separately the indices of both zero and nonzero val‐
762 ues in the mask.
763 where returns associated values from a data PDL, rather than indices
765 into the mask PDL.
766 whichND returns N-D indices into a multidimensional PDL.
768 perldl> $x = sequence(10); p $x
770 [0 1 2 3 4 5 6 7 8 9]
771 perldl> $indx = which($x>6); p $indx
772 [7 8 9]
773 which_both
775 Signature: (mask(n); int [o] inds(m); int [o]notinds(q))
777 Returns indices of zero and nonzero values in a mask PDL
779 ($i, $c_i) = which_both($mask);
781 This works just as which, but the complement of $i will be in $c_i.
783 perldl> $x = sequence(10); p $x
785 [0 1 2 3 4 5 6 7 8 9]
786 perldl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
787 [5 6 7 8 9]
788 [0 1 2 3 4]
789 where
791 Use a mask to select values from one or more data PDLs
793 "where" accepts one or more data piddles and a mask piddle. It returns
795 a list of output piddles, corresponding to the input data piddles.
796 Each output piddle is a 1-dimensional list of values in its correspond‐
797 ing data piddle. The values are drawn from locations where the mask is
798 nonzero.
799 The output PDLs are still connected to the original data PDLs, for the
801 purpose of dataflow.
802 "where" combines the functionality of which and index into a single
804 operation.
805 BUGS:
807 There is no "whereND", and probably should be. While "where" works OK
809 for most N-dimensional cases, it does not thread properly over (for
810 example) the (N+1)th dimension in data that is compared to an N-dimen‐
811 sional mask.
812 $i = $x->where($x+5 > 0); # $i contains those elements of $x
814 # where mask ($x+5 > 0) is 1
815 $i .= -5; # Set those elements (of $x) to -5. Together, these
816 # commands clamp $x to a maximum of -5.
817 It is also possible to use the same mask for several piddles with the
819 same call:
820 ($i,$j,$k) = where($x,$y,$z, $x+5>0);
822 Note: $i is always 1-D, even if $x is >1-D.
824 WARNING: The first argument (the values) and the second argument (the
826 mask) currently have to have the exact same dimensions (or horrible
827 things happen). You *cannot* thread over a smaller mask, for example.
828 whichND
830 Return the coordinates of non-zero values in a mask.
832 WhichND returns the N-dimensional coordinates of each nonzero value in
834 a mask PDL with any number of dimensions.
835 For historical reasons the return value is different in list and scalar
837 context. In scalar context, you get back a PDL containing coordinates
838 suitable for use in indexND or range; in list context, the coordinates
839 are broken out into separate PDLs.
840 $coords = whichND($mask);
842 returns a PDL containing the coordinates of the elements that are non-
844 zero in $mask, suitable for use in indexND. The 0th dimension contains
845 the full coordinate listing of each point; the 1st dimension lists all
846 the points. For example, if $mask has rank 4 and 100 matching ele‐
847 ments, then $coords has dimension 4x100.
848 @coords=whichND($mask);
850 returns a perl list of piddles containing the coordinates of the ele‐
852 ments that are non-zero in $mask. Each element corresponds to a par‐
853 ticular index dimension. For example, if $mask has rank 4 and 100
854 matching elements, then @coords has 4 elements, each of which is a pdl
855 of size 100.
856 SEE ALSO:
858 which finds coordinates of nonzero values in a 1-D mask.
860 where extracts values from a data PDL that are associated with nonzero
862 values in a mask PDL.
863 perldl> $a=sequence(10,10,3,4)
865 perldl> ($x, $y, $z, $w)=whichND($a == 203); p $x, $y, $z, $w
866 [3] [0] [2] [0]
867 perldl> print $a->at(list(cat($x,$y,$z,$w)))
868 203
869 setops
871 Implements simple set operations like union and intersection
873 Usage: $set = setops($a, <OPERATOR>, $b);
875 The operator can be "OR", "XOR" or "AND". This is then applied to $a
877 viewed as a set and $b viewed as a set. The functioning is as follows:
878 "OR"
880 The resulting vector will contain the elements that are either in
881 $a or in $b or both. This is the union in set operation terms
882 "XOR"
884 The resulting vector will contain the elements that are either in
885 $a or $b, but not in both. This is
886 Union($a, $b) - Intersection($a, $b)
888 in set operation terms.
890 "AND"
892 The resulting vector will contain the intersection of $a and $b, so
893 the elements that are in both $a and $b. Note that for convenience
894 this operation is also aliased to intersect
895 It should be emphasized that these routines are used when one or both
897 of the sets $a, $b are hard to calculate or that you get from a sepa‐
898 rate subroutine.
899 Finally IDL users might be familiar with Craig Markwardt's
901 "cmset_op.pro" routine which has inspired this routine although it was
902 written independently However the present routine has a few less
903 options (but see the exampels)
904 You will very often use these functions on an index vector, so that is
906 what we will show here. We will in fact something slightly silly. First
907 we will find all squares that are also cubes below 10000.
908 Create a sequence vector:
910 perldl> $x = sequence(10000)
912 Find all odd and even elements:
914 perldl> ($even, $odd) = which_both( ($x % 2) == 0)
916 Find all squares
918 perldl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
920 Find all cubes (being careful with roundoff error!)
922 perldl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
924 Then find all squares that are cubes:
926 perldl> $both = setops($squares, 'AND', $cubes)
928 And print these (assumes that "PDL::NiceSlice" is loaded!)
930 perldl> p $x($both)
932 [0 1 64 729 4096]
933 Then find all numbers that are either cubes or squares, but not both:
935 perldl> $cube_xor_square = setops($squares, 'XOR', $cubes)
937 perldl> p $cube_xor_square->nelem()
939 112
940 So there are a total of 112 of these!
942 Finally find all odd squares:
944 perldl> $odd_squares = setops($squares, 'AND', $odd)
946 Another common occurance is to want to get all objects that are in $a
948 and in the complement of $b. But it is almost always best to create the
949 complement explicitly since the universe that both are taken from is
950 not known. Thus use which_both if possible to keep track of comple‐
951 ments.
952 If this is impossible the best approach is to make a temporary:
954 This creates an index vector the size of the universe of the sets and
956 set all elements in $b to 0
957 perldl> $tmp = ones($n_universe); $tmp($b)=0;
959 This then finds the complement of $b
961 perldl> $C_b = which($tmp == 1);
963 and this does the final selection:
965 perldl> $set = setops($a, 'AND', $C_b)
967 intersect
969 Calculate the intersection of two piddles
971 Usage: $set = intersect($a, $b);
973 This routine is merely a simple interface to setops. See that for more
975 information
976 Find all numbers less that 100 that are of the form 2*y and 3*x
978 perldl> $x=sequence(100)
980 perldl> $factor2 = which( ($x % 2) == 0)
982 perldl> $factor3 = which( ($x % 3) == 0)
984 perldl> $ii=intersect($factor2, $factor3)
986 perldl> p $x($ii)
988 [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
990
992 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu). Contribu‐
993 tions by Christian Soeller (c.soeller@auckland.ac.nz), Karl Glazebrook
994 (kgb@aaoepp.aao.gov.au), Craig DeForest (deforest@boulder.swri.edu) and
995 Jarle Brinchmann (jarle@astro.up.pt) All rights reserved. There is no
996 warranty. You are allowed to redistribute this software / documentation
997 under certain conditions. For details, see the file COPYING in the PDL
998 distribution. If this file is separated from the PDL distribution, the
999 copyright notice should be included in the file.
1000perl v5.8.8 2006-12-02 Primitive(3)