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())
21 Inner product over one dimension
23 c = sum_i a_i * b_i
25 If "a() * b()" contains only bad data, "c()" is set bad. Otherwise
27 "c()" will have its bad flag cleared, as it will not contain any bad
28 values.
29 outer
31 Signature: (a(n); b(m); [o]c(n,m))
32 outer product over one dimension
34 Naturally, it is possible to achieve the effects of outer product
36 simply by threading over the ""*"" operator but this function is
37 provided for convenience.
38 outer does handle bad values. It will set the bad-value flag of all
40 output piddles if the flag is set for any of the input piddles.
41 x
43 Signature: (a(i,x), b(z,i),[o]c(x,z))
44 Matrix multiplication
46 PDL overloads the "x" operator (normally the repeat operator) for
48 matrix multiplication. The number of columns (size of the 0 dimension)
49 in the left-hand argument must normally equal the number of rows (size
50 of the 1 dimension) in the right-hand argument.
51 Row vectors are represented as (N x 1) two-dimensional PDLs, or you may
53 be sloppy and use a one-dimensional PDL. Column vectors are
54 represented as (1 x N) two-dimensional PDLs.
55 Threading occurs in the usual way, but as both the 0 and 1 dimension
57 (if present) are included in the operation, you must be sure that you
58 don't try to thread over either of those dims.
59 EXAMPLES
61 Here are some simple ways to define vectors and matrices:
63 perldl> $r = pdl(1,2); # A row vector
65 perldl> $c = pdl([[3],[4]]); # A column vector
66 perldl> $c = pdl(3,4)->(*1); # A column vector, using NiceSlice
67 perldl> $m = pdl([[1,2],[3,4]]); # A 2x2 matrix
68 Now that we have a few objects prepared, here is how to matrix-multiply
70 them:
71 perldl> print $r x $m # row x matrix = row
73 [
74 [ 7 10]
75 ]
76 perldl> print $m x $r # matrix x row = ERROR
78 PDL: Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1
79 perldl> print $m x $c # matrix x column = column
81 [
82 [ 5]
83 [11]
84 ]
85 perldl> print $m x 2 # Trivial case: scalar mult.
87 [
88 [2 4]
89 [6 8]
90 ]
91 perldl> print $r x $c # row x column = scalar
93 [
94 [11]
95 ]
96 perldl> print $c x $r # column x row = matrix
98 [
99 [3 6]
100 [4 8]
101 ]
102 INTERNALS
104 The mechanics of the multiplication are carried out by the matmult
106 method.
107 matmult
109 Signature: (a(i,x),b(z,i),[o]c(x,z))
110 Matrix multiplication
112 We peruse the inner product to define matrix multiplication via a
114 threaded inner product.
115 For usage, see x, a description of the overloaded 'x' operator
117 innerwt
119 Signature: (a(n); b(n); c(n); [o]d())
120 Weighted (i.e. triple) inner product
122 d = sum_i a(i) b(i) c(i)
124 innerwt does handle bad values. It will set the bad-value flag of all
126 output piddles if the flag is set for any of the input piddles.
127 inner2
129 Signature: (a(n); b(n,m); c(m); [o]d())
130 Inner product of two vectors and a matrix
132 d = sum_ij a(i) b(i,j) c(j)
134 Note that you should probably not thread over "a" and "c" since that
136 would be very wasteful. Instead, you should use a temporary for "b*c".
137 inner2 does handle bad values. It will set the bad-value flag of all
139 output piddles if the flag is set for any of the input piddles.
140 inner2d
142 Signature: (a(n,m); b(n,m); [o]c())
143 Inner product over 2 dimensions.
145 Equivalent to
147 $c = inner($a->clump(2), $b->clump(2))
149 inner2d does handle bad values. It will set the bad-value flag of all
151 output piddles if the flag is set for any of the input piddles.
152 inner2t
154 Signature: (a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k)))
155 Efficient Triple matrix product "a*b*c"
157 Efficiency comes from by using the temporary "tmp". This operation only
159 scales as "N**3" whereas threading using inner2 would scale as "N**4".
160 The reason for having this routine is that you do not need to have the
162 same thread-dimensions for "tmp" as for the other arguments, which in
163 case of large numbers of matrices makes this much more memory-
164 efficient.
165 It is hoped that things like this could be taken care of as a kind of
167 closures at some point.
168 inner2t does handle bad values. It will set the bad-value flag of all
170 output piddles if the flag is set for any of the input piddles.
171 crossp
173 Signature: (a(tri=3); b(tri); [o] c(tri))
174 Cross product of two 3D vectors
176 After
178 $c = crossp $a, $b
180 the inner product "$c*$a" and "$c*$b" will be zero, i.e. $c is
182 orthogonal to $a and $b
183 crossp does not process bad values. It will set the bad-value flag of
185 all output piddles if the flag is set for any of the input piddles.
186 norm
188 Signature: (vec(n); [o] norm(n))
189 Normalises a vector to unit Euclidean length
191 norm does handle bad values. It will set the bad-value flag of all
193 output piddles if the flag is set for any of the input piddles.
194 indadd
196 Signature: (a(); int ind(); [o] sum(m))
197 Threaded Index Add: Add "a" to the "ind" element of "sum", i.e:
199 sum(ind) += a
201 Simple Example:
203 $a = 2;
205 $ind = 3;
206 $sum = zeroes(10);
207 indadd($a,$ind, $sum);
208 print $sum
209 #Result: ( 2 added to element 3 of $sum)
210 # [0 0 0 2 0 0 0 0 0 0]
211 Threaded Example:
213 $a = pdl( 1,2,3);
215 $ind = pdl( 1,4,6);
216 $sum = zeroes(10);
217 indadd($a,$ind, $sum);
218 print $sum."\n";
219 #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
220 # [0 1 0 0 2 0 3 0 0 0]
221 The routine barfs if any of the indices are bad.
223 conv1d
225 Signature: (a(m); kern(p); [o]b(m); int reflect)
226 1D convolution along first dimension
228 The m-th element of the discrete convolution of an input piddle $a of
230 size $M, and a kernel piddle $kern of size $P, is calculated as
231 n = ($P-1)/2
233 ====
234 \
235 ($a conv1d $kern)[m] = > $a_ext[m - n] * $kern[n]
236 /
237 ====
238 n = -($P-1)/2
239 where $a_ext is either the periodic (or reflected) extension of $a so
241 it is equal to $a on " 0..$M-1 " and equal to the corresponding
242 periodic/reflected image of $a outside that range.
243 $con = conv1d sequence(10), pdl(-1,0,1);
245 $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
247 By default, periodic boundary conditions are assumed (i.e. wrap
249 around). Alternatively, you can request reflective boundary conditions
250 using the "Boundary" option:
251 {Boundary => 'reflect'} # case in 'reflect' doesn't matter
253 The convolution is performed along the first dimension. To apply it
255 across another dimension use the slicing routines, e.g.
256 $b = $a->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
258 This function is useful for threaded filtering of 1D signals.
260 Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv
262 conv1d does not process bad values. It will set the bad-value flag of
264 all output piddles if the flag is set for any of the input piddles.
265 in
267 Signature: (a(); b(n); [o] c())
268 test if a is in the set of values b
270 $goodmsk = $labels->in($goodlabels);
272 print pdl(4,3,1)->in(pdl(2,3,3));
273 [0 1 0]
274 "in" is akin to the is an element of of set theory. In priciple, PDL
276 threading could be used to achieve its functionality by using a
277 construct like
278 $msk = ($labels->dummy(0) == $goodlabels)->orover;
280 However, "in" doesn't create a (potentially large) intermediate and is
282 generally faster.
283 in does not process bad values. It will set the bad-value flag of all
285 output piddles if the flag is set for any of the input piddles.
286 uniq
288 return all unique elements of a piddle
289 The unique elements are returned in ascending order.
291 print pdl(2,2,2,4,0,-1,6,6)->uniq;
293 [-1 0 2 4 6]
294 Note: The returned pdl is 1D; any structure of the input piddle is
296 lost.
297 See uniqind if you need the indices of the unique elements rather than
299 the values.
300 Bad values are not considered unique by uniq and are ignored.
302 $a=sequence(10);
304 $a=$a->setbadif($a%3);
305 print $a->uniq;
306 [0 3 6 9]
307 uniqind
309 return the indices of all unique elements of a piddle The order is in
310 the order of the values to be consistent with uniq
311 print pdl(2,2,2,4,0,-1,6,6)->uniqind;
313 [5, 4, 1, 3, 6]
314 Note: The returned pdl is 1D; any structure of the input piddle is
316 lost.
317 See uniq if you want the unique values instead of the indices.
319 Bad values are not considered unique by uniqind and are ignored.
321 uniqvec
323 return all unique vectors out of a collection
324 The unique vectors are returned in lexicographically sorted ascending
326 order. The 0th dimension of the input PDL is treated as a dimensional
327 index within each vector, and the 1st and any higher dimensions are
328 taken to run across vectors. The return value is always 2D; any
329 structure of the input PDL (beyond using the 0th dimension for vector
330 index) is lost.
331 See also uniq for a uniqe list of scalars; and qsortvec for sorting a
333 list of vectors lexicographcally.
334 If a vector contains all bad values, it is ignored as in uniq. If some
336 of the values are good, it is treated as a normal vector. For example,
337 [1 2 BAD] and [BAD 2 3] could be returned, but [BAD BAD BAD] could not.
338 hclip
340 Signature: (a(); b(); [o] c())
341 clip (threshold) $a by $b ($b is upper bound)
343 hclip does handle bad values. It will set the bad-value flag of all
345 output piddles if the flag is set for any of the input piddles.
346 lclip
348 Signature: (a(); b(); [o] c())
349 clip (threshold) $a by $b ($b is lower bound)
351 lclip does handle bad values. It will set the bad-value flag of all
353 output piddles if the flag is set for any of the input piddles.
354 clip
356 Clip (threshold) a piddle by (optional) upper or lower bounds.
357 $b = $a->clip(0,3);
359 $c = $a->clip(undef, $x);
360 clip handles bad values since it is just a wrapper around hclip and
362 lclip.
363 wtstat
365 Signature: (a(n); wt(n); avg(); [o]b(); int deg)
366 Weighted statistical moment of given degree
368 This calculates a weighted statistic over the vector "a". The formula
370 is
371 b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
373 Bad values are ignored in any calculation; $b will only have its bad
375 flag set if the output contains any bad data.
376 statsover
378 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())
379 Calculate useful statistics over a dimension of a piddle
381 ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($piddle, $weights);
383 This utility function calculates various useful quantities of a piddle.
385 These are:
386 · the mean:
388 MEAN = sum (x)/ N
390 with "N" being the number of elements in x
392 · RMS deviation from the mean:
394 RMS = sqrt(sum( (x-mean(x))^2 )/N)
396 (also known as the root-mean-square deviation, or the square root of
398 the variance)
399 · the median
401 The median is the 50th percentile data value. Median is found by
403 medover, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
404 · the minimum
406 · the maximum
408 · the absolute deviation:
410 ADEV = sqrt(sum( abs(x-mean(x)) )/N)
412 (This is also called the standard deviation)
414 · the population RMS deviation from the mean:
416 PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
418 The population deviation is the best-estimate of the deviation of
420 the population from which a sample is drawn.
421 This operator is a projection operator so the calculation will take
423 place over the final dimension. Thus if the input is N-dimensional each
424 returned value will be N-1 dimensional, to calculate the statistics for
425 the entire piddle either use "clump(-1)" directly on the piddle or call
426 "stats".
427 Bad values are simply ignored in the calculation, effectively reducing
429 the sample size. If all data are bad then the output data are marked
430 bad.
431 stats
433 Calculates useful statistics on a piddle
434 ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($piddle,[$weights]);
436 This utility calculates all the most useful quantities in one call. It
438 works the same way as "statsover", except that the quantities are
439 calculated considering the entire input PDL as a single sample, rather
440 than as a collection of rows. See "statsover" for definitions of the
441 returned quantities.
442 Bad values are handled; if all input values are bad, then all of the
444 output values are flagged bad.
445 histogram
447 Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
448 Calculates a histogram for given stepsize and minimum.
450 $h = histogram($data, $step, $min, $numbins);
452 $hist = zeroes $numbins; # Put histogram in existing piddle.
453 histogram($data, $hist, $step, $min, $numbins);
454 The histogram will contain $numbins bins starting from $min, each $step
456 wide. The value in each bin is the number of values in $data that lie
457 within the bin limits.
458 Data below the lower limit is put in the first bin, and data above the
460 upper limit is put in the last bin.
461 The output is reset in a different threadloop so that you can take a
463 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
464 Use hist instead for a high-level interface.
466 perldl> p histogram(pdl(1,1,2),1,0,3)
468 [0 2 1]
469 histogram does handle bad values. It will set the bad-value flag of
471 all output piddles if the flag is set for any of the input piddles.
472 whistogram
474 Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
475 Calculates a histogram from weighted data for given stepsize and
477 minimum.
478 $h = whistogram($data, $weights, $step, $min, $numbins);
480 $hist = zeroes $numbins; # Put histogram in existing piddle.
481 whistogram($data, $weights, $hist, $step, $min, $numbins);
482 The histogram will contain $numbins bins starting from $min, each $step
484 wide. The value in each bin is the sum of the values in $weights that
485 correspond to values in $data that lie within the bin limits.
486 Data below the lower limit is put in the first bin, and data above the
488 upper limit is put in the last bin.
489 The output is reset in a different threadloop so that you can take a
491 histogram of "$a(10,12)" into "$b(15)" and get the result you want.
492 perldl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
494 [0 0.2 0.5 0]
495 whistogram does handle bad values. It will set the bad-value flag of
497 all output piddles if the flag is set for any of the input piddles.
498 histogram2d
500 Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
501 double stepb; double minb; int mbsize => mb;)
502 Calculates a 2d histogram.
504 $h = histogram2d($datax, $datay,
506 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
507 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
508 histogram2d($datax, $datay, $hist,
509 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
510 The histogram will contain $nbinx x $nbiny bins, with the lower limits
512 of the first one at "($minx, $miny)", and with bin size "($stepx,
513 $stepy)". The value in each bin is the number of values in $datax and
514 $datay that lie within the bin limits.
515 Data below the lower limit is put in the first bin, and data above the
517 upper limit is put in the last bin.
518 perldl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
520 [
521 [0 0 0]
522 [0 2 2]
523 [0 1 0]
524 ]
525 histogram2d does handle bad values. It will set the bad-value flag of
527 all output piddles if the flag is set for any of the input piddles.
528 whistogram2d
530 Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
531 double stepb; double minb; int mbsize => mb;)
532 Calculates a 2d histogram from weighted data.
534 $h = whistogram2d($datax, $datay, $weights,
536 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
537 $hist = zeroes $nbinx, $nbiny; # Put histogram in existing piddle.
538 whistogram2d($datax, $datay, $weights, $hist,
539 $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
540 The histogram will contain $nbinx x $nbiny bins, with the lower limits
542 of the first one at "($minx, $miny)", and with bin size "($stepx,
543 $stepy)". The value in each bin is the sum of the values in $weights
544 that correspond to values in $datax and $datay that lie within the bin
545 limits.
546 Data below the lower limit is put in the first bin, and data above the
548 upper limit is put in the last bin.
549 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)
551 [
552 [ 0 0 0]
553 [ 0 0.5 0.9]
554 [ 0 0.1 0]
555 ]
556 whistogram2d does handle bad values. It will set the bad-value flag of
558 all output piddles if the flag is set for any of the input piddles.
559 fibonacci
561 Signature: ([o]x(n))
562 Constructor - a vector with Fibonacci's sequence
564 fibonacci does not process bad values. It will set the bad-value flag
566 of all output piddles if the flag is set for any of the input piddles.
567 append
569 Signature: (a(n); b(m); [o] c(mn))
570 append two or more piddles by concatenating along their first
572 dimensions
573 $a = ones(2,4,7);
575 $b = sequence 5;
576 $c = $a->append($b); # size of $c is now (7,4,7) (a jumbo-piddle ;)
577 "append" appends two piddles along their first dims. Rest of the
579 dimensions must be compatible in the threading sense. Resulting size of
580 first dim is the sum of the sizes of the first dims of the two argument
581 piddles - ie "n + m".
582 append does not process bad values. It will set the bad-value flag of
584 all output piddles if the flag is set for any of the input piddles.
585 glue
587 $c = $a->glue(<dim>,$b,...)
588 Glue two or more PDLs together along an arbitrary dimension (N-D
590 append).
591 Sticks $a, $b, and all following arguments together along the specified
593 dimension. All other dimensions must be compatible in the threading
594 sense.
595 Glue is permissive, in the sense that every PDL is treated as having an
597 infinite number of trivial dimensions of order 1 -- so "$a-"glue(3,$b)>
598 works, even if $a and $b are only one dimensional.
599 If one of the PDLs has no elements, it is ignored. Likewise, if one of
601 them is actually the undefined value, it is treated as if it had no
602 elements.
603 If the first parameter is a defined perl scalar rather than a pdl, then
605 it is taken as a dimension along which to glue everything else, so you
606 can say "$cube = PDL::glue(3,@image_list);" if you like.
607 "glue" is implemented in pdl, using a combination of xchg and append.
609 It should probably be updated (one day) to a pure PP function.
610 axisvalues
612 Signature: ([o,nc]a(n))
613 Internal routine
615 "axisvalues" is the internal primitive that implements axisvals and
617 alters its argument.
618 axisvalues does not process bad values. It will set the bad-value flag
620 of all output piddles if the flag is set for any of the input piddles.
621 random
623 Constructor which returns piddle of random numbers
624 $a = random([type], $nx, $ny, $nz,...);
626 $a = random $b;
627 etc (see zeroes).
629 This is the uniform distribution between 0 and 1 (assumedly excluding 1
631 itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
632 specify dimensions, types or give a template.
633 You can use the perl function srand to seed the random generator. For
635 further details consult Perl's srand documentation.
636 randsym
638 Constructor which returns piddle of random numbers
639 $a = randsym([type], $nx, $ny, $nz,...);
641 $a = randsym $b;
642 etc (see zeroes).
644 This is the uniform distribution between 0 and 1 (excluding both 0 and
646 1, cf random). The arguments are the same as "zeroes" (q.v.) - i.e. one
647 can specify dimensions, types or give a template.
648 You can use the perl function srand to seed the random generator. For
650 further details consult Perl's srand documentation.
651 grandom
653 Constructor which returns piddle of Gaussian random numbers
654 $a = grandom([type], $nx, $ny, $nz,...);
656 $a = grandom $b;
657 etc (see zeroes).
659 This is generated using the math library routine "ndtri".
661 Mean = 0, Stddev = 1
663 You can use the perl function srand to seed the random generator. For
665 further details consult Perl's srand documentation.
666 vsearch
668 Signature: (i(); x(n); int [o]ip())
669 routine for searching 1D values i.e. step-function interpolation.
671 $inds = vsearch($vals, $xs);
673 Returns for each value of $vals the index of the least larger member of
675 $xs (which need to be in increasing order). If the value is larger than
676 any member of $xs, the index to the last element of $xs is returned.
677 This function is useful e.g. when you have a list of probabilities for
679 events and want to generate indices to events:
680 $a = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
682 $b = random 20;
683 $c = vsearch($b, $a); # Now, $c will have the appropriate distr.
684 It is possible to use the cumusumover function to obtain cumulative
686 probabilities from absolute probabilities.
687 needs major (?) work to handles bad values
689 interpolate
691 Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
692 routine for 1D linear interpolation
694 ( $yi, $err ) = interpolate($xi, $x, $y)
696 Given a set of points "($x,$y)", use linear interpolation to find the
698 values $yi at a set of points $xi.
699 "interpolate" uses a binary search to find the suspects, er...,
701 interpolation indices and therefore abscissas (ie $x) have to be
702 strictly ordered (increasing or decreasing). For interpolation at lots
703 of closely spaced abscissas an approach that uses the last index found
704 as a start for the next search can be faster (compare Numerical Recipes
705 "hunt" routine). Feel free to implement that on top of the binary
706 search if you like. For out of bounds values it just does a linear
707 extrapolation and sets the corresponding element of $err to 1, which is
708 otherwise 0.
709 See also interpol, which uses the same routine, differing only in the
711 handling of extrapolation - an error message is printed rather than
712 returning an error piddle.
713 needs major (?) work to handles bad values
715 interpol
717 Signature: (xi(); x(n); y(n); [o] yi())
718 routine for 1D linear interpolation
720 $yi = interpol($xi, $x, $y)
722 "interpol" uses the same search method as interpolate, hence $x must be
724 strictly ordered (either increasing or decreasing). The difference
725 occurs in the handling of out-of-bounds values; here an error message
726 is printed.
727 interpND
729 Interpolate values from an N-D piddle, with switchable method
730 $source = 10*xvals(10,10) + yvals(10,10);
732 $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
733 print $source->interpND( $index );
734 InterpND acts like indexND, collapsing $index by lookup into $source;
736 but it does interpolation rather than direct sampling. The
737 interpolation method and boundary condition are switchable via an
738 options hash.
739 By default, linear or sample interpolation is used, with constant value
741 outside the boundaries of the source pdl. No dataflow occurs, because
742 in general the output is computed rather than indexed.
743 All the interpolation methods treat the pixels as value-centered, so
745 the "sample" method will return $a->(0) for coordinate values on the
746 set [-0.5,0.5), and all methods will return $a->(1) for a coordinate
747 value of exactly 1.
748 Recognized options:
750 method
752 Values can be:
753 · 0, s, sample, Sample (default for integer source types)
755 The nearest value is taken. Pixels are regarded as centered on
757 their respective integer coordinates (no offset from the linear
758 case).
759 · 1, l, linear, Linear (default for floating point source types)
761 The values are N-linearly interpolated from an N-dimensional cube
763 of size 2.
764 · 3, c, cube, cubic, Cubic
766 The values are interpolated using a local cubic fit to the data.
768 The fit is constrained to match the original data and its
769 derivative at the data points. The second derivative of the fit
770 is not continuous at the data points. Multidimensional datasets
771 are interpolated by the successive-collapse method.
772 (Note that the constraint on the first derivative causes a small
774 amount of ringing around sudden features such as step functions).
775 · f, fft, fourier, Fourier
777 The source is Fourier transformed, and the interpolated values
779 are explicitly calculated from the coefficients. The boundary
780 condition option is ignored -- periodic boundaries are imposed.
781 If you pass in the option "fft", and it is a list (ARRAY) ref,
783 then it is a stash for the magnitude and phase of the source FFT.
784 If the list has two elements then they are taken as already
785 computed; otherwise they are calculated and put in the stash.
786 b, bound, boundary, Boundary
788 This option is passed unmodified into indexND, which is used as the
789 indexing engine for the interpolation. Some current allowed values
790 are 'extend', 'periodic', 'truncate', and 'mirror' (default is
791 'truncate').
792 bad
794 contains the fill value used for 'truncate' boundary. (default 0)
795 fft
797 An array ref whose associated list is used to stash the FFT of the
798 source data, for the FFT method.
799 one2nd
801 Converts a one dimensional index piddle to a set of ND coordinates
802 @coords=one2nd($a, $indices)
804 returns an array of piddles containing the ND indexes corresponding to
806 the one dimensional list indices. The indices are assumed to correspond
807 to array $a clumped using "clump(-1)". This routine is used in whichND,
808 but is useful on its own occasionally.
809 perldl> $a=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$a->clump(-1)
811 perldl> $maxind=maximum_ind($c); p $maxind;
812 6
813 perldl> print one2nd($a, maximum_ind($c))
814 0 1 1
815 perldl> p $a->at(0,1,1)
816 3
817 which
819 Signature: (mask(n); int [o] inds(m))
820 Returns indices of non-zero values from a 1-D PDL
822 $i = which($mask);
824 returns a pdl with indices for all those elements that are nonzero in
826 the mask. Note that the returned indices will be 1D. If you feed in a
827 multidimensional mask, it will be flattened before the indices are
828 calculated. See also whichND for multidimensional masks.
829 If you want to index into the original mask or a similar piddle with
831 output from "which", remember to flatten it before calling index:
832 $data = random 5, 5;
834 $idx = which $data > 0.5; # $idx is now 1D
835 $bigsum = $data->flat->index($idx)->sum; # flatten before indexing
836 Compare also where for similar functionality.
838 SEE ALSO:
840 which_both returns separately the indices of both zero and nonzero
842 values in the mask.
843 where returns associated values from a data PDL, rather than indices
845 into the mask PDL.
846 whichND returns N-D indices into a multidimensional PDL.
848 perldl> $x = sequence(10); p $x
850 [0 1 2 3 4 5 6 7 8 9]
851 perldl> $indx = which($x>6); p $indx
852 [7 8 9]
853 which does handle bad values. It will set the bad-value flag of all
855 output piddles if the flag is set for any of the input piddles.
856 which_both
858 Signature: (mask(n); int [o] inds(m); int [o]notinds(q))
859 Returns indices of zero and nonzero values in a mask PDL
861 ($i, $c_i) = which_both($mask);
863 This works just as which, but the complement of $i will be in $c_i.
865 perldl> $x = sequence(10); p $x
867 [0 1 2 3 4 5 6 7 8 9]
868 perldl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
869 [5 6 7 8 9]
870 [0 1 2 3 4]
871 which_both does handle bad values. It will set the bad-value flag of
873 all output piddles if the flag is set for any of the input piddles.
874 where
876 Use a mask to select values from one or more data PDLs
877 "where" accepts one or more data piddles and a mask piddle. It returns
879 a list of output piddles, corresponding to the input data piddles.
880 Each output piddle is a 1-dimensional list of values in its
881 corresponding data piddle. The values are drawn from locations where
882 the mask is nonzero.
883 The output PDLs are still connected to the original data PDLs, for the
885 purpose of dataflow.
886 "where" combines the functionality of which and index into a single
888 operation.
889 BUGS:
891 There is no "whereND", and probably should be. While "where" works OK
893 for most N-dimensional cases, it does not thread properly over (for
894 example) the (N+1)th dimension in data that is compared to an
895 N-dimensional mask.
896 $i = $x->where($x+5 > 0); # $i contains those elements of $x
898 # where mask ($x+5 > 0) is 1
899 $i .= -5; # Set those elements (of $x) to -5. Together, these
900 # commands clamp $x to a maximum of -5.
901 It is also possible to use the same mask for several piddles with the
903 same call:
904 ($i,$j,$k) = where($x,$y,$z, $x+5>0);
906 Note: $i is always 1-D, even if $x is >1-D.
908 WARNING: The first argument (the values) and the second argument (the
910 mask) currently have to have the exact same dimensions (or horrible
911 things happen). You *cannot* thread over a smaller mask, for example.
912 whichND
914 Return the coordinates of non-zero values in a mask.
915 WhichND returns the N-dimensional coordinates of each nonzero value in
917 a mask PDL with any number of dimensions.
918 For historical reasons the return value is different in list and scalar
920 context. In scalar context, you get back a PDL containing coordinates
921 suitable for use in indexND or range; in list context, the coordinates
922 are broken out into separate PDLs.
923 $coords = whichND($mask);
925 returns a PDL containing the coordinates of the elements that are non-
927 zero in $mask, suitable for use in indexND. The 0th dimension contains
928 the full coordinate listing of each point; the 1st dimension lists all
929 the points. For example, if $mask has rank 4 and 100 matching
930 elements, then $coords has dimension 4x100.
931 @coords=whichND($mask);
933 returns a perl list of piddles containing the coordinates of the
935 elements that are non-zero in $mask. Each element corresponds to a
936 particular index dimension. For example, if $mask has rank 4 and 100
937 matching elements, then @coords has 4 elements, each of which is a pdl
938 of size 100.
939 SEE ALSO:
941 which finds coordinates of nonzero values in a 1-D mask.
943 where extracts values from a data PDL that are associated with nonzero
945 values in a mask PDL.
946 perldl> $a=sequence(10,10,3,4)
948 perldl> ($x, $y, $z, $w)=whichND($a == 203); p $x, $y, $z, $w
949 [3] [0] [2] [0]
950 perldl> print $a->at(list(cat($x,$y,$z,$w)))
951 203
952 setops
954 Implements simple set operations like union and intersection
955 Usage: $set = setops($a, <OPERATOR>, $b);
957 The operator can be "OR", "XOR" or "AND". This is then applied to $a
959 viewed as a set and $b viewed as a set. Set theory says that a set may
960 not have two or more identical elements, but setops takes care of this
961 for you, so "$a=pdl(1,1,2)" is OK. The functioning is as follows:
962 "OR"
964 The resulting vector will contain the elements that are either in
965 $a or in $b or both. This is the union in set operation terms
966 "XOR"
968 The resulting vector will contain the elements that are either in
969 $a or $b, but not in both. This is
970 Union($a, $b) - Intersection($a, $b)
972 in set operation terms.
974 "AND"
976 The resulting vector will contain the intersection of $a and $b, so
977 the elements that are in both $a and $b. Note that for convenience
978 this operation is also aliased to intersect
979 It should be emphasized that these routines are used when one or both
981 of the sets $a, $b are hard to calculate or that you get from a
982 separate subroutine.
983 Finally IDL users might be familiar with Craig Markwardt's
985 "cmset_op.pro" routine which has inspired this routine although it was
986 written independently However the present routine has a few less
987 options (but see the exampels)
988 You will very often use these functions on an index vector, so that is
990 what we will show here. We will in fact something slightly silly. First
991 we will find all squares that are also cubes below 10000.
992 Create a sequence vector:
994 perldl> $x = sequence(10000)
996 Find all odd and even elements:
998 perldl> ($even, $odd) = which_both( ($x % 2) == 0)
1000 Find all squares
1002 perldl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
1004 Find all cubes (being careful with roundoff error!)
1006 perldl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
1008 Then find all squares that are cubes:
1010 perldl> $both = setops($squares, 'AND', $cubes)
1012 And print these (assumes that "PDL::NiceSlice" is loaded!)
1014 perldl> p $x($both)
1016 [0 1 64 729 4096]
1017 Then find all numbers that are either cubes or squares, but not both:
1019 perldl> $cube_xor_square = setops($squares, 'XOR', $cubes)
1021 perldl> p $cube_xor_square->nelem()
1023 112
1024 So there are a total of 112 of these!
1026 Finally find all odd squares:
1028 perldl> $odd_squares = setops($squares, 'AND', $odd)
1030 Another common occurance is to want to get all objects that are in $a
1032 and in the complement of $b. But it is almost always best to create the
1033 complement explicitly since the universe that both are taken from is
1034 not known. Thus use which_both if possible to keep track of
1035 complements.
1036 If this is impossible the best approach is to make a temporary:
1038 This creates an index vector the size of the universe of the sets and
1040 set all elements in $b to 0
1041 perldl> $tmp = ones($n_universe); $tmp($b)=0;
1043 This then finds the complement of $b
1045 perldl> $C_b = which($tmp == 1);
1047 and this does the final selection:
1049 perldl> $set = setops($a, 'AND', $C_b)
1051 intersect
1053 Calculate the intersection of two piddles
1054 Usage: $set = intersect($a, $b);
1056 This routine is merely a simple interface to setops. See that for more
1058 information
1059 Find all numbers less that 100 that are of the form 2*y and 3*x
1061 perldl> $x=sequence(100)
1063 perldl> $factor2 = which( ($x % 2) == 0)
1065 perldl> $factor3 = which( ($x % 3) == 0)
1067 perldl> $ii=intersect($factor2, $factor3)
1069 perldl> p $x($ii)
1071 [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
1073
1075 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
1076 Contributions by Christian Soeller (c.soeller@auckland.ac.nz), Karl
1077 Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
1078 (deforest@boulder.swri.edu) and Jarle Brinchmann (jarle@astro.up.pt)
1079 All rights reserved. There is no warranty. You are allowed to
1080 redistribute this software / documentation under certain conditions.
1081 For details, see the file COPYING in the PDL distribution. If this file
1082 is separated from the PDL distribution, the copyright notice should be
1083 included in the file.
1084perl v5.12.3 2011-03-31 Primitive(3)