1Primitive(3)          User Contributed Perl Documentation         Primitive(3)
2
3
4

NAME

6       PDL::Primitive - primitive operations for pdl
7

DESCRIPTION

9       This module provides some primitive and useful functions defined using
10       PDL::PP and able to use the new indexing tricks.
11
12       See PDL::Indexing for how to use indices creatively.  For explanation
13       of the signature format, see PDL::PP.
14

SYNOPSIS

16        # Pulls in PDL::Primitive, among other modules.
17        use PDL;
18
19        # Only pull in PDL::Primitive:
20        use PDL::Primitive;
21

FUNCTIONS

23   inner
24         Signature: (a(n); b(n); [o]c())
25
26       Inner product over one dimension
27
28        c = sum_i a_i * b_i
29
30       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
34   outer
35         Signature: (a(n); b(m); [o]c(n,m))
36
37       outer product over one dimension
38
39       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
43       outer processes bad values.  It will set the bad-value flag of all
44       output piddles if the flag is set for any of the input piddles.
45
46   x
47        Signature: (a(i,z), b(x,i),[o]c(x,z))
48
49       Matrix multiplication
50
51       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
56       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
60       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
64       EXAMPLES
65
66       Here are some simple ways to define vectors and matrices:
67
68        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
73       Now that we have a few objects prepared, here is how to matrix-multiply
74       them:
75
76        pdl> print $r x $m                 # row x matrix = row
77        [
78         [ 7 10]
79        ]
80
81        pdl> print $m x $r                 # matrix x row = ERROR
82        PDL: Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1
83
84        pdl> print $m x $c                 # matrix x column = column
85        [
86         [ 5]
87         [11]
88        ]
89
90        pdl> print $m x 2                  # Trivial case: scalar mult.
91        [
92         [2 4]
93         [6 8]
94        ]
95
96        pdl> print $r x $c                 # row x column = scalar
97        [
98         [11]
99        ]
100
101        pdl> print $c x $r                 # column x row = matrix
102        [
103         [3 6]
104         [4 8]
105        ]
106
107       INTERNALS
108
109       The mechanics of the multiplication are carried out by the matmult
110       method.
111
112   matmult
113         Signature: (a(t,h); b(w,t); [o]c(w,h))
114
115       Matrix multiplication
116
117       Notionally, matrix multiplication $a x $b is equivalent to the
118       threading expression
119
120           $a->dummy(1)->inner($b->xchg(0,1)->dummy(2),$c);
121
122       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
126       For usage, see x, a description of the overloaded 'x' operator
127
128       matmult ignores the bad-value flag of the input piddles.  It will set
129       the bad-value flag of all output piddles if the flag is set for any of
130       the input piddles.
131
132   innerwt
133         Signature: (a(n); b(n); c(n); [o]d())
134
135       Weighted (i.e. triple) inner product
136
137        d = sum_i a(i) b(i) c(i)
138
139       innerwt processes bad values.  It will set the bad-value flag of all
140       output piddles if the flag is set for any of the input piddles.
141
142   inner2
143         Signature: (a(n); b(n,m); c(m); [o]d())
144
145       Inner product of two vectors and a matrix
146
147        d = sum_ij a(i) b(i,j) c(j)
148
149       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
152       inner2 processes bad values.  It will set the bad-value flag of all
153       output piddles if the flag is set for any of the input piddles.
154
155   inner2d
156         Signature: (a(n,m); b(n,m); [o]c())
157
158       Inner product over 2 dimensions.
159
160       Equivalent to
161
162        $c = inner($a->clump(2), $b->clump(2))
163
164       inner2d processes bad values.  It will set the bad-value flag of all
165       output piddles if the flag is set for any of the input piddles.
166
167   inner2t
168         Signature: (a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k)))
169
170       Efficient Triple matrix product "a*b*c"
171
172       Efficiency comes from by using the temporary "tmp". This operation only
173       scales as "N**3" whereas threading using inner2 would scale as "N**4".
174
175       The reason for having this routine is that you do not need to have the
176       same thread-dimensions for "tmp" as for the other arguments, which in
177       case of large numbers of matrices makes this much more memory-
178       efficient.
179
180       It is hoped that things like this could be taken care of as a kind of
181       closures at some point.
182
183       inner2t processes bad values.  It will set the bad-value flag of all
184       output piddles if the flag is set for any of the input piddles.
185
186   crossp
187         Signature: (a(tri=3); b(tri); [o] c(tri))
188
189       Cross product of two 3D vectors
190
191       After
192
193        $c = crossp $a, $b
194
195       the inner product "$c*$a" and "$c*$b" will be zero, i.e. $c is
196       orthogonal to $a and $b
197
198       crossp does not process bad values.  It will set the bad-value flag of
199       all output piddles if the flag is set for any of the input piddles.
200
201   norm
202         Signature: (vec(n); [o] norm(n))
203
204       Normalises a vector to unit Euclidean length
205
206       norm processes bad values.  It will set the bad-value flag of all
207       output piddles if the flag is set for any of the input piddles.
208
209   indadd
210         Signature: (a(); indx ind(); [o] sum(m))
211
212       Threaded Index Add: Add "a" to the "ind" element of "sum", i.e:
213
214        sum(ind) += a
215
216       Simple Example:
217
218         $a = 2;
219         $ind = 3;
220         $sum = zeroes(10);
221         indadd($a,$ind, $sum);
222         print $sum
223         #Result: ( 2 added to element 3 of $sum)
224         # [0 0 0 2 0 0 0 0 0 0]
225
226       Threaded Example:
227
228         $a = pdl( 1,2,3);
229         $ind = pdl( 1,4,6);
230         $sum = zeroes(10);
231         indadd($a,$ind, $sum);
232         print $sum."\n";
233         #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
234         # [0 1 0 0 2 0 3 0 0 0]
235
236       The routine barfs if any of the indices are bad.
237
238   conv1d
239         Signature: (a(m); kern(p); [o]b(m); int reflect)
240
241       1D convolution along first dimension
242
243       The m-th element of the discrete convolution of an input piddle $a of
244       size $M, and a kernel piddle $kern of size $P, is calculated as
245
246                                     n = ($P-1)/2
247                                     ====
248                                     \
249         ($a conv1d $kern)[m]   =     >      $a_ext[m - n] * $kern[n]
250                                     /
251                                     ====
252                                     n = -($P-1)/2
253
254       where $a_ext is either the periodic (or reflected) extension of $a so
255       it is equal to $a on " 0..$M-1 " and equal to the corresponding
256       periodic/reflected image of $a outside that range.
257
258         $con = conv1d sequence(10), pdl(-1,0,1);
259
260         $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
261
262       By default, periodic boundary conditions are assumed (i.e. wrap
263       around).  Alternatively, you can request reflective boundary conditions
264       using the "Boundary" option:
265
266         {Boundary => 'reflect'} # case in 'reflect' doesn't matter
267
268       The convolution is performed along the first dimension. To apply it
269       across another dimension use the slicing routines, e.g.
270
271         $b = $a->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
272
273       This function is useful for threaded filtering of 1D signals.
274
275       Compare also conv2d, convolve, fftconvolve, fftwconv, rfftwconv
276
277       WARNING: "conv1d" processes bad values in its inputs as the numeric
278       value of "$pdl->badvalue" so it is not recommended for processing pdls
279       with bad values in them unless special care is taken.
280
281       conv1d ignores the bad-value flag of the input piddles.  It will set
282       the bad-value flag of all output piddles if the flag is set for any of
283       the input piddles.
284
285   in
286         Signature: (a(); b(n); [o] c())
287
288       test if a is in the set of values b
289
290          $goodmsk = $labels->in($goodlabels);
291          print pdl(3,1,4,6,2)->in(pdl(2,3,3));
292         [1 0 0 0 1]
293
294       "in" is akin to the is an element of of set theory. In principle, PDL
295       threading could be used to achieve its functionality by using a
296       construct like
297
298          $msk = ($labels->dummy(0) == $goodlabels)->orover;
299
300       However, "in" doesn't create a (potentially large) intermediate and is
301       generally faster.
302
303       in does not process bad values.  It will set the bad-value flag of all
304       output piddles if the flag is set for any of the input piddles.
305
306   uniq
307       return all unique elements of a piddle
308
309       The unique elements are returned in ascending order.
310
311         PDL> p pdl(2,2,2,4,0,-1,6,6)->uniq
312         [-1 0 2 4 6]     # 0 is returned 2nd (sorted order)
313
314         PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniq
315         [-1 2 4 6 nan]   # NaN value is returned at end
316
317       Note: The returned pdl is 1D; any structure of the input piddle is
318       lost.  "NaN" values are never compare equal to any other values, even
319       themselves.  As a result, they are always unique. "uniq" returns the
320       NaN values at the end of the result piddle.  This follows the Matlab
321       usage.
322
323       See uniqind if you need the indices of the unique elements rather than
324       the values.
325
326       Bad values are not considered unique by uniq and are ignored.
327
328        $a=sequence(10);
329        $a=$a->setbadif($a%3);
330        print $a->uniq;
331        [0 3 6 9]
332
333   uniqind
334       Return the indices of all unique elements of a piddle The order is in
335       the order of the values to be consistent with uniq. "NaN" values never
336       compare equal with any other value and so are always unique.  This
337       follows the Matlab usage.
338
339         PDL> p pdl(2,2,2,4,0,-1,6,6)->uniqind
340         [5 4 1 3 6]     # the 0 at index 4 is returned 2nd, but...
341
342         PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniqind
343         [5 1 3 6 4]     # ...the NaN at index 4 is returned at end
344
345       Note: The returned pdl is 1D; any structure of the input piddle is
346       lost.
347
348       See uniq if you want the unique values instead of the indices.
349
350       Bad values are not considered unique by uniqind and are ignored.
351
352   uniqvec
353       Return all unique vectors out of a collection
354
355         NOTE: If any vectors in the input piddle have NaN values
356         they are returned at the end of the non-NaN ones.  This is
357         because, by definition, NaN values never compare equal with
358         any other value.
359
360         NOTE: The current implementation does not sort the vectors
361         containing NaN values.
362
363       The unique vectors are returned in lexicographically sorted ascending
364       order. The 0th dimension of the input PDL is treated as a dimensional
365       index within each vector, and the 1st and any higher dimensions are
366       taken to run across vectors. The return value is always 2D; any
367       structure of the input PDL (beyond using the 0th dimension for vector
368       index) is lost.
369
370       See also uniq for a unique list of scalars; and qsortvec for sorting a
371       list of vectors lexicographcally.
372
373       If a vector contains all bad values, it is ignored as in uniq.  If some
374       of the values are good, it is treated as a normal vector. For example,
375       [1 2 BAD] and [BAD 2 3] could be returned, but [BAD BAD BAD] could not.
376       Vectors containing BAD values will be returned after any non-NaN and
377       non-BAD containing vectors, followed by the NaN vectors.
378
379   hclip
380         Signature: (a(); b(); [o] c())
381
382       clip (threshold) $a by $b ($b is upper bound)
383
384       hclip processes bad values.  It will set the bad-value flag of all
385       output piddles if the flag is set for any of the input piddles.
386
387   lclip
388         Signature: (a(); b(); [o] c())
389
390       clip (threshold) $a by $b ($b is lower bound)
391
392       lclip processes bad values.  It will set the bad-value flag of all
393       output piddles if the flag is set for any of the input piddles.
394
395   clip
396       Clip (threshold) a piddle by (optional) upper or lower bounds.
397
398        $b = $a->clip(0,3);
399        $c = $a->clip(undef, $x);
400
401       clip handles bad values since it is just a wrapper around hclip and
402       lclip.
403
404   clip
405         Signature: (a(); l(); h(); [o] c())
406
407       info not available
408
409       clip processes bad values.  It will set the bad-value flag of all
410       output piddles if the flag is set for any of the input piddles.
411
412   wtstat
413         Signature: (a(n); wt(n); avg(); [o]b(); int deg)
414
415       Weighted statistical moment of given degree
416
417       This calculates a weighted statistic over the vector "a".  The formula
418       is
419
420        b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
421
422       Bad values are ignored in any calculation; $b will only have its bad
423       flag set if the output contains any bad data.
424
425   statsover
426         Signature: (a(n); w(n); float+ [o]avg(); float+ [o]prms(); int+ [o]median(); int+ [o]min(); int+ [o]max(); float+ [o]adev(); float+ [o]rms())
427
428       Calculate useful statistics over a dimension of a piddle
429
430         ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($piddle, $weights);
431
432       This utility function calculates various useful quantities of a piddle.
433       These are:
434
435       ·  the mean:
436
437            MEAN = sum (x)/ N
438
439          with "N" being the number of elements in x
440
441       ·  the population RMS deviation from the mean:
442
443            PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
444
445          The population deviation is the best-estimate of the deviation of
446          the population from which a sample is drawn.
447
448       ·  the median
449
450          The median is the 50th percentile data value.  Median is found by
451          medover, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
452
453       ·  the minimum
454
455       ·  the maximum
456
457       ·  the average absolute deviation:
458
459            AADEV = sum( abs(x-mean(x)) )/N
460
461       ·  RMS deviation from the mean:
462
463            RMS = sqrt(sum( (x-mean(x))^2 )/N)
464
465          (also known as the root-mean-square deviation, or the square root of
466          the variance)
467
468       This operator is a projection operator so the calculation will take
469       place over the final dimension. Thus if the input is N-dimensional each
470       returned value will be N-1 dimensional, to calculate the statistics for
471       the entire piddle either use "clump(-1)" directly on the piddle or call
472       "stats".
473
474       Bad values are simply ignored in the calculation, effectively reducing
475       the sample size.  If all data are bad then the output data are marked
476       bad.
477
478   stats
479       Calculates useful statistics on a piddle
480
481        ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($piddle,[$weights]);
482
483       This utility calculates all the most useful quantities in one call.  It
484       works the same way as "statsover", except that the quantities are
485       calculated considering the entire input PDL as a single sample, rather
486       than as a collection of rows. See "statsover" for definitions of the
487       returned quantities.
488
489       Bad values are handled; if all input values are bad, then all of the
490       output values are flagged bad.
491
492   histogram
493         Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
494
495       Calculates a histogram for given stepsize and minimum.
496
497        $h = histogram($data, $step, $min, $numbins);
498        $hist = zeroes $numbins;  # Put histogram in existing piddle.
499        histogram($data, $hist, $step, $min, $numbins);
500
501       The histogram will contain $numbins bins starting from $min, each $step
502       wide. The value in each bin is the number of values in $data that lie
503       within the bin limits.
504
505       Data below the lower limit is put in the first bin, and data above the
506       upper limit is put in the last bin.
507
508       The output is reset in a different threadloop so that you can take a
509       histogram of "$a(10,12)" into "$b(15)" and get the result you want.
510
511       For a higher-level interface, see hist.
512
513        pdl> p histogram(pdl(1,1,2),1,0,3)
514        [0 2 1]
515
516       histogram processes bad values.  It will set the bad-value flag of all
517       output piddles if the flag is set for any of the input piddles.
518
519   whistogram
520         Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
521
522       Calculates a histogram from weighted data for given stepsize and
523       minimum.
524
525        $h = whistogram($data, $weights, $step, $min, $numbins);
526        $hist = zeroes $numbins;  # Put histogram in existing piddle.
527        whistogram($data, $weights, $hist, $step, $min, $numbins);
528
529       The histogram will contain $numbins bins starting from $min, each $step
530       wide. The value in each bin is the sum of the values in $weights that
531       correspond to values in $data that lie within the bin limits.
532
533       Data below the lower limit is put in the first bin, and data above the
534       upper limit is put in the last bin.
535
536       The output is reset in a different threadloop so that you can take a
537       histogram of "$a(10,12)" into "$b(15)" and get the result you want.
538
539        pdl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
540        [0 0.2 0.5 0]
541
542       whistogram processes bad values.  It will set the bad-value flag of all
543       output piddles if the flag is set for any of the input piddles.
544
545   histogram2d
546         Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
547                            double stepb; double minb; int mbsize => mb;)
548
549       Calculates a 2d histogram.
550
551        $h = histogram2d($datax, $datay, $stepx, $minx,
552              $nbinx, $stepy, $miny, $nbiny);
553        $hist = zeroes $nbinx, $nbiny;  # Put histogram in existing piddle.
554        histogram2d($datax, $datay, $hist, $stepx, $minx,
555              $nbinx, $stepy, $miny, $nbiny);
556
557       The histogram will contain $nbinx x $nbiny bins, with the lower limits
558       of the first one at "($minx, $miny)", and with bin size "($stepx,
559       $stepy)".  The value in each bin is the number of values in $datax and
560       $datay that lie within the bin limits.
561
562       Data below the lower limit is put in the first bin, and data above the
563       upper limit is put in the last bin.
564
565        pdl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
566        [
567         [0 0 0]
568         [0 2 2]
569         [0 1 0]
570        ]
571
572       histogram2d processes bad values.  It will set the bad-value flag of
573       all output piddles if the flag is set for any of the input piddles.
574
575   whistogram2d
576         Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
577                            double stepb; double minb; int mbsize => mb;)
578
579       Calculates a 2d histogram from weighted data.
580
581        $h = whistogram2d($datax, $datay, $weights,
582              $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
583        $hist = zeroes $nbinx, $nbiny;  # Put histogram in existing piddle.
584        whistogram2d($datax, $datay, $weights, $hist,
585              $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
586
587       The histogram will contain $nbinx x $nbiny bins, with the lower limits
588       of the first one at "($minx, $miny)", and with bin size "($stepx,
589       $stepy)".  The value in each bin is the sum of the values in $weights
590       that correspond to values in $datax and $datay that lie within the bin
591       limits.
592
593       Data below the lower limit is put in the first bin, and data above the
594       upper limit is put in the last bin.
595
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)
597        [
598         [  0   0   0]
599         [  0 0.5 0.9]
600         [  0 0.1   0]
601        ]
602
603       whistogram2d processes bad values.  It will set the bad-value flag of
604       all output piddles if the flag is set for any of the input piddles.
605
606   fibonacci
607         Signature: ([o]x(n))
608
609       Constructor - a vector with Fibonacci's sequence
610
611       fibonacci does not process bad values.  It will set the bad-value flag
612       of all output piddles if the flag is set for any of the input piddles.
613
614   append
615         Signature: (a(n); b(m); [o] c(mn))
616
617       append two piddles by concatenating along their first dimensions
618
619        $a = ones(2,4,7);
620        $b = sequence 5;
621        $c = $a->append($b);  # size of $c is now (7,4,7) (a jumbo-piddle ;)
622
623       "append" appends two piddles along their first dimensions. The rest of
624       the dimensions must be compatible in the threading sense. The resulting
625       size of the first dimension is the sum of the sizes of the first
626       dimensions of the two argument piddles - i.e. "n + m".
627
628       Similar functions include glue (below), which can append more than two
629       piddles along an arbitary dimension, and cat, which can append more
630       than two piddles that all have the same sized dimensions.
631
632       append does not process bad values.  It will set the bad-value flag of
633       all output piddles if the flag is set for any of the input piddles.
634
635   glue
636         $c = $a->glue(<dim>,$b,...)
637
638       Glue two or more PDLs together along an arbitrary dimension (N-D
639       append).
640
641       Sticks $a, $b, and all following arguments together along the specified
642       dimension.  All other dimensions must be compatible in the threading
643       sense.
644
645       Glue is permissive, in the sense that every PDL is treated as having an
646       infinite number of trivial dimensions of order 1 -- so "$a->glue(3,$b)"
647       works, even if $a and $b are only one dimensional.
648
649       If one of the PDLs has no elements, it is ignored.  Likewise, if one of
650       them is actually the undefined value, it is treated as if it had no
651       elements.
652
653       If the first parameter is a defined perl scalar rather than a pdl, then
654       it is taken as a dimension along which to glue everything else, so you
655       can say "$cube = PDL::glue(3,@image_list);" if you like.
656
657       "glue" is implemented in pdl, using a combination of xchg and append.
658       It should probably be updated (one day) to a pure PP function.
659
660       Similar functions include append (above), which appends only two
661       piddles along their first dimension, and cat, which can append more
662       than two piddles that all have the same sized dimensions.
663
664   axisvalues
665         Signature: ([o,nc]a(n))
666
667       Internal routine
668
669       "axisvalues" is the internal primitive that implements axisvals and
670       alters its argument.
671
672       axisvalues does not process bad values.  It will set the bad-value flag
673       of all output piddles if the flag is set for any of the input piddles.
674
675   random
676       Constructor which returns piddle of random numbers
677
678        $a = random([type], $nx, $ny, $nz,...);
679        $a = random $b;
680
681       etc (see zeroes).
682
683       This is the uniform distribution between 0 and 1 (assumedly excluding 1
684       itself). The arguments are the same as "zeroes" (q.v.) - i.e. one can
685       specify dimensions, types or give a template.
686
687       You can use the perl function srand to seed the random generator. For
688       further details consult Perl's  srand documentation.
689
690   randsym
691       Constructor which returns piddle of random numbers
692
693        $a = randsym([type], $nx, $ny, $nz,...);
694        $a = randsym $b;
695
696       etc (see zeroes).
697
698       This is the uniform distribution between 0 and 1 (excluding both 0 and
699       1, cf random). The arguments are the same as "zeroes" (q.v.) - i.e. one
700       can specify dimensions, types or give a template.
701
702       You can use the perl function srand to seed the random generator. For
703       further details consult Perl's  srand documentation.
704
705   grandom
706       Constructor which returns piddle of Gaussian random numbers
707
708        $a = grandom([type], $nx, $ny, $nz,...);
709        $a = grandom $b;
710
711       etc (see zeroes).
712
713       This is generated using the math library routine "ndtri".
714
715       Mean = 0, Stddev = 1
716
717       You can use the perl function srand to seed the random generator. For
718       further details consult Perl's  srand documentation.
719
720   vsearch
721         Signature: ( vals(); xs(n); [o] indx(); [\%options] )
722
723       Efficiently search for values in a sorted piddle, returning indices.
724
725         $idx = vsearch( $vals, $x, [\%options] );
726         vsearch( $vals, $x, $idx, [\%options ] );
727
728       vsearch performs a binary search in the ordered piddle $x, for the
729       values from $vals piddle, returning indices into $x.  What is a
730       "match", and the meaning of the returned indices, are determined by the
731       options.
732
733       The "mode" option indicates which method of searching to use, and may
734       be one of:
735
736       "sample"
737           invoke vsearch_sample, returning indices appropriate for sampling
738           within a distribution.
739
740       "insert_leftmost"
741           invoke vsearch_insert_leftmost, returning the left-most possible
742           insertion point which still leaves the piddle sorted.
743
744       "insert_rightmost"
745           invoke vsearch_insert_rightmost, returning the right-most possible
746           insertion point which still leaves the piddle sorted.
747
748       "insert_match"
749           invoke vsearch_match, returning the index of a matching element,
750           else -(insertion point + 1)
751
752       "insert_bin_inclusive"
753           invoke vsearch_bin_inclusive, returning an index appropriate for
754           binning on a grid where the left bin edges are inclusive of the
755           bin. See below for further explanation of the bin.
756
757       "insert_bin_exclusive"
758           invoke vsearch_bin_exclusive, returning an index appropriate for
759           binning on a grid where the left bin edges are exclusive of the
760           bin. See below for further explanation of the bin.
761
762       The default value of "mode" is "sample".
763
764   vsearch_sample
765         Signature: (vals(); x(n); indx [o]idx())
766
767       Search for values in a sorted array, return index appropriate for
768       sampling from a distribution
769
770         $idx = vsearch_sample($vals, $x);
771
772       $x must be sorted, but may be in decreasing or increasing order.
773
774       vsearch_sample returns an index I for each value V of $vals appropriate
775       for sampling $vals
776
777       I has the following properties:
778
779       ·   if $x is sorted in increasing order
780
781                     V <= x[0]  : I = 0
782             x[0]  < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
783             x[-1] < V          : I = $x->nelem -1
784
785       ·   if $x is sorted in decreasing order
786
787                      V > x[0]  : I = 0
788             x[0]  >= V > x[-1] : I s.t. x[I] >= V > x[I+1]
789             x[-1] >= V         : I = $x->nelem - 1
790
791       If all elements of $x are equal, I = $x->nelem - 1.
792
793       If $x contains duplicated elements, I is the index of the leftmost (by
794       position in array) duplicate if V matches.
795
796       This function is useful e.g. when you have a list of probabilities for
797       events and want to generate indices to events:
798
799        $a = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
800        $b = random 20;
801        $c = vsearch_sample($b, $a); # Now, $c will have the appropriate distr.
802
803       It is possible to use the cumusumover function to obtain cumulative
804       probabilities from absolute probabilities.
805
806       needs major (?) work to handles bad values
807
808   vsearch_insert_leftmost
809         Signature: (vals(); x(n); indx [o]idx())
810
811       Determine the insertion point for values in a sorted array, inserting
812       before duplicates.
813
814         $idx = vsearch_insert_leftmost($vals, $x);
815
816       $x must be sorted, but may be in decreasing or increasing order.
817
818       vsearch_insert_leftmost returns an index I for each value V of $vals
819       equal to the leftmost position (by index in array) within $x that V may
820       be inserted and still maintain the order in $x.
821
822       Insertion at index I involves shifting elements I and higher of $x to
823       the right by one and setting the now empty element at index I to V.
824
825       I has the following properties:
826
827       ·   if $x is sorted in increasing order
828
829                     V <= x[0]  : I = 0
830             x[0]  < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
831             x[-1] < V          : I = $x->nelem
832
833       ·   if $x is sorted in decreasing order
834
835                      V >  x[0]  : I = -1
836             x[0]  >= V >= x[-1] : I s.t. x[I] >= V > x[I+1]
837             x[-1] >= V          : I = $x->nelem -1
838
839       If all elements of $x are equal,
840
841           i = 0
842
843       If $x contains duplicated elements, I is the index of the leftmost (by
844       index in array) duplicate if V matches.
845
846       needs major (?) work to handles bad values
847
848   vsearch_insert_rightmost
849         Signature: (vals(); x(n); indx [o]idx())
850
851       Determine the insertion point for values in a sorted array, inserting
852       after duplicates.
853
854         $idx = vsearch_insert_rightmost($vals, $x);
855
856       $x must be sorted, but may be in decreasing or increasing order.
857
858       vsearch_insert_rightmost returns an index I for each value V of $vals
859       equal to the rightmost position (by index in array) within $x that V
860       may be inserted and still maintain the order in $x.
861
862       Insertion at index I involves shifting elements I and higher of $x to
863       the right by one and setting the now empty element at index I to V.
864
865       I has the following properties:
866
867       ·   if $x is sorted in increasing order
868
869                      V < x[0]  : I = 0
870             x[0]  <= V < x[-1] : I s.t. x[I-1] <= V < x[I]
871             x[-1] <= V         : I = $x->nelem
872
873       ·   if $x is sorted in decreasing order
874
875                     V >= x[0]  : I = -1
876             x[0]  > V >= x[-1] : I s.t. x[I] >= V > x[I+1]
877             x[-1] > V          : I = $x->nelem -1
878
879       If all elements of $x are equal,
880
881           i = $x->nelem - 1
882
883       If $x contains duplicated elements, I is the index of the leftmost (by
884       index in array) duplicate if V matches.
885
886       needs major (?) work to handles bad values
887
888   vsearch_match
889         Signature: (vals(); x(n); indx [o]idx())
890
891       Match values against a sorted array.
892
893         $idx = vsearch_match($vals, $x);
894
895       $x must be sorted, but may be in decreasing or increasing order.
896
897       vsearch_match returns an index I for each value V of $vals.  If V
898       matches an element in $x, I is the index of that element, otherwise it
899       is -( insertion_point + 1 ), where insertion_point is an index in $x
900       where V may be inserted while maintaining the order in $x.  If $x has
901       duplicated values, I may refer to any of them.
902
903       needs major (?) work to handles bad values
904
905   vsearch_bin_inclusive
906         Signature: (vals(); x(n); indx [o]idx())
907
908       Determine the index for values in a sorted array of bins, lower bound
909       inclusive.
910
911         $idx = vsearch_bin_inclusive($vals, $x);
912
913       $x must be sorted, but may be in decreasing or increasing order.
914
915       $x represents the edges of contiguous bins, with the first and last
916       elements representing the outer edges of the outer bins, and the inner
917       elements the shared bin edges.
918
919       The lower bound of a bin is inclusive to the bin, its outer bound is
920       exclusive to it.  vsearch_bin_inclusive returns an index I for each
921       value V of $vals
922
923       I has the following properties:
924
925       ·   if $x is sorted in increasing order
926
927                      V < x[0]  : I = -1
928             x[0]  <= V < x[-1] : I s.t. x[I] <= V < x[I+1]
929             x[-1] <= V         : I = $x->nelem - 1
930
931       ·   if $x is sorted in decreasing order
932
933                      V >= x[0]  : I = 0
934             x[0]  >  V >= x[-1] : I s.t. x[I+1] > V >= x[I]
935             x[-1] >  V          : I = $x->nelem
936
937       If all elements of $x are equal,
938
939           i = $x->nelem - 1
940
941       If $x contains duplicated elements, I is the index of the righmost (by
942       index in array) duplicate if V matches.
943
944       needs major (?) work to handles bad values
945
946   vsearch_bin_exclusive
947         Signature: (vals(); x(n); indx [o]idx())
948
949       Determine the index for values in a sorted array of bins, lower bound
950       exclusive.
951
952         $idx = vsearch_bin_exclusive($vals, $x);
953
954       $x must be sorted, but may be in decreasing or increasing order.
955
956       $x represents the edges of contiguous bins, with the first and last
957       elements representing the outer edges of the outer bins, and the inner
958       elements the shared bin edges.
959
960       The lower bound of a bin is exclusive to the bin, its upper bound is
961       inclusive to it.  vsearch_bin_exclusive returns an index I for each
962       value V of $vals.
963
964       I has the following properties:
965
966       ·   if $x is sorted in increasing order
967
968                      V <= x[0]  : I = -1
969             x[0]  <  V <= x[-1] : I s.t. x[I] < V <= x[I+1]
970             x[-1] <  V          : I = $x->nelem - 1
971
972       ·   if $x is sorted in decreasing order
973
974                      V >  x[0]  : I = 0
975             x[0]  >= V >  x[-1] : I s.t. x[I-1] >= V > x[I]
976             x[-1] >= V          : I = $x->nelem
977
978       If all elements of $x are equal,
979
980           i = $x->nelem - 1
981
982       If $x contains duplicated elements, I is the index of the righmost (by
983       index in array) duplicate if V matches.
984
985       needs major (?) work to handles bad values
986
987   interpolate
988         Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
989
990       routine for 1D linear interpolation
991
992        ( $yi, $err ) = interpolate($xi, $x, $y)
993
994       Given a set of points "($x,$y)", use linear interpolation to find the
995       values $yi at a set of points $xi.
996
997       "interpolate" uses a binary search to find the suspects, er...,
998       interpolation indices and therefore abscissas (ie $x) have to be
999       strictly ordered (increasing or decreasing).  For interpolation at lots
1000       of closely spaced abscissas an approach that uses the last index found
1001       as a start for the next search can be faster (compare Numerical Recipes
1002       "hunt" routine). Feel free to implement that on top of the binary
1003       search if you like. For out of bounds values it just does a linear
1004       extrapolation and sets the corresponding element of $err to 1, which is
1005       otherwise 0.
1006
1007       See also interpol, which uses the same routine, differing only in the
1008       handling of extrapolation - an error message is printed rather than
1009       returning an error piddle.
1010
1011       needs major (?) work to handles bad values
1012
1013   interpol
1014        Signature: (xi(); x(n); y(n); [o] yi())
1015
1016       routine for 1D linear interpolation
1017
1018        $yi = interpol($xi, $x, $y)
1019
1020       "interpol" uses the same search method as interpolate, hence $x must be
1021       strictly ordered (either increasing or decreasing).  The difference
1022       occurs in the handling of out-of-bounds values; here an error message
1023       is printed.
1024
1025   interpND
1026       Interpolate values from an N-D piddle, with switchable method
1027
1028         $source = 10*xvals(10,10) + yvals(10,10);
1029         $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
1030         print $source->interpND( $index );
1031
1032       InterpND acts like indexND, collapsing $index by lookup into $source;
1033       but it does interpolation rather than direct sampling.  The
1034       interpolation method and boundary condition are switchable via an
1035       options hash.
1036
1037       By default, linear or sample interpolation is used, with constant value
1038       outside the boundaries of the source pdl.  No dataflow occurs, because
1039       in general the output is computed rather than indexed.
1040
1041       All the interpolation methods treat the pixels as value-centered, so
1042       the "sample" method will return "$a->(0)" for coordinate values on the
1043       set [-0.5,0.5), and all methods will return "$a->(1)" for a coordinate
1044       value of exactly 1.
1045
1046       Recognized options:
1047
1048       method
1049          Values can be:
1050
1051          ·  0, s, sample, Sample (default for integer source types)
1052
1053             The nearest value is taken. Pixels are regarded as centered on
1054             their respective integer coordinates (no offset from the linear
1055             case).
1056
1057          ·  1, l, linear, Linear (default for floating point source types)
1058
1059             The values are N-linearly interpolated from an N-dimensional cube
1060             of size 2.
1061
1062          ·  3, c, cube, cubic, Cubic
1063
1064             The values are interpolated using a local cubic fit to the data.
1065             The fit is constrained to match the original data and its
1066             derivative at the data points.  The second derivative of the fit
1067             is not continuous at the data points.  Multidimensional datasets
1068             are interpolated by the successive-collapse method.
1069
1070             (Note that the constraint on the first derivative causes a small
1071             amount of ringing around sudden features such as step functions).
1072
1073          ·  f, fft, fourier, Fourier
1074
1075             The source is Fourier transformed, and the interpolated values
1076             are explicitly calculated from the coefficients.  The boundary
1077             condition option is ignored -- periodic boundaries are imposed.
1078
1079             If you pass in the option "fft", and it is a list (ARRAY) ref,
1080             then it is a stash for the magnitude and phase of the source FFT.
1081             If the list has two elements then they are taken as already
1082             computed; otherwise they are calculated and put in the stash.
1083
1084       b, bound, boundary, Boundary
1085          This option is passed unmodified into indexND, which is used as the
1086          indexing engine for the interpolation.  Some current allowed values
1087          are 'extend', 'periodic', 'truncate', and 'mirror' (default is
1088          'truncate').
1089
1090       bad
1091          contains the fill value used for 'truncate' boundary.  (default 0)
1092
1093       fft
1094          An array ref whose associated list is used to stash the FFT of the
1095          source data, for the FFT method.
1096
1097   one2nd
1098       Converts a one dimensional index piddle to a set of ND coordinates
1099
1100        @coords=one2nd($a, $indices)
1101
1102       returns an array of piddles containing the ND indexes corresponding to
1103       the one dimensional list indices. The indices are assumed to correspond
1104       to array $a clumped using "clump(-1)". This routine is used in the old
1105       vector form of whichND, but is useful on its own occasionally.
1106
1107       Returned piddles have the indx datatype.  $indices can have values
1108       larger than "$a->nelem" but negative values in $indices will not give
1109       the answer you expect.
1110
1111        pdl> $a=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$a->clump(-1)
1112        pdl> $maxind=maximum_ind($c); p $maxind;
1113        6
1114        pdl> print one2nd($a, maximum_ind($c))
1115        0 1 1
1116        pdl> p $a->at(0,1,1)
1117        3
1118
1119   which
1120         Signature: (mask(n); indx [o] inds(m))
1121
1122       Returns indices of non-zero values from a 1-D PDL
1123
1124        $i = which($mask);
1125
1126       returns a pdl with indices for all those elements that are nonzero in
1127       the mask. Note that the returned indices will be 1D. If you feed in a
1128       multidimensional mask, it will be flattened before the indices are
1129       calculated.  See also whichND for multidimensional masks.
1130
1131       If you want to index into the original mask or a similar piddle with
1132       output from "which", remember to flatten it before calling index:
1133
1134         $data = random 5, 5;
1135         $idx = which $data > 0.5; # $idx is now 1D
1136         $bigsum = $data->flat->index($idx)->sum;  # flatten before indexing
1137
1138       Compare also where for similar functionality.
1139
1140       SEE ALSO:
1141
1142       which_both returns separately the indices of both zero and nonzero
1143       values in the mask.
1144
1145       where returns associated values from a data PDL, rather than indices
1146       into the mask PDL.
1147
1148       whichND returns N-D indices into a multidimensional PDL.
1149
1150        pdl> $x = sequence(10); p $x
1151        [0 1 2 3 4 5 6 7 8 9]
1152        pdl> $indx = which($x>6); p $indx
1153        [7 8 9]
1154
1155       which processes bad values.  It will set the bad-value flag of all
1156       output piddles if the flag is set for any of the input piddles.
1157
1158   which_both
1159         Signature: (mask(n); indx [o] inds(m); indx [o]notinds(q))
1160
1161       Returns indices of zero and nonzero values in a mask PDL
1162
1163        ($i, $c_i) = which_both($mask);
1164
1165       This works just as which, but the complement of $i will be in $c_i.
1166
1167        pdl> $x = sequence(10); p $x
1168        [0 1 2 3 4 5 6 7 8 9]
1169        pdl> ($small, $big) = which_both ($x >= 5); p "$small\n $big"
1170        [5 6 7 8 9]
1171        [0 1 2 3 4]
1172
1173       which_both processes bad values.  It will set the bad-value flag of all
1174       output piddles if the flag is set for any of the input piddles.
1175
1176   where
1177       Use a mask to select values from one or more data PDLs
1178
1179       "where" accepts one or more data piddles and a mask piddle.  It returns
1180       a list of output piddles, corresponding to the input data piddles.
1181       Each output piddle is a 1-dimensional list of values in its
1182       corresponding data piddle. The values are drawn from locations where
1183       the mask is nonzero.
1184
1185       The output PDLs are still connected to the original data PDLs, for the
1186       purpose of dataflow.
1187
1188       "where" combines the functionality of which and index into a single
1189       operation.
1190
1191       BUGS:
1192
1193       While "where" works OK for most N-dimensional cases, it does not thread
1194       properly over (for example) the (N+1)th dimension in data that is
1195       compared to an N-dimensional mask.  Use "whereND" for that.
1196
1197        $i = $x->where($x+5 > 0); # $i contains those elements of $x
1198                                  # where mask ($x+5 > 0) is 1
1199        $i .= -5;  # Set those elements (of $x) to -5. Together, these
1200                   # commands clamp $x to a maximum of -5.
1201
1202       It is also possible to use the same mask for several piddles with the
1203       same call:
1204
1205        ($i,$j,$k) = where($x,$y,$z, $x+5>0);
1206
1207       Note: $i is always 1-D, even if $x is >1-D.
1208
1209       WARNING: The first argument (the values) and the second argument (the
1210       mask) currently have to have the exact same dimensions (or horrible
1211       things happen). You *cannot* thread over a smaller mask, for example.
1212
1213   whereND
1214       "where" with support for ND masks and threading
1215
1216       "whereND" accepts one or more data piddles and a mask piddle.  It
1217       returns a list of output piddles, corresponding to the input data
1218       piddles.  The values are drawn from locations where the mask is
1219       nonzero.
1220
1221       "whereND" differs from "where" in that the mask dimensionality is
1222       preserved which allows for proper threading of the selection operation
1223       over higher dimensions.
1224
1225       As with "where" the output PDLs are still connected to the original
1226       data PDLs, for the purpose of dataflow.
1227
1228         $sdata = whereND $data, $mask
1229         ($s1, $s2, ..., $sn) = whereND $d1, $d2, ..., $dn, $mask
1230
1231         where
1232
1233           $data is M dimensional
1234           $mask is N < M dimensional
1235           dims($data) 1..N == dims($mask) 1..N
1236           with threading over N+1 to M dimensions
1237
1238         $data   = sequence(4,3,2);   # example data array
1239         $mask4  = (random(4)>0.5);   # example 1-D mask array, has $n4 true values
1240         $mask43 = (random(4,3)>0.5); # example 2-D mask array, has $n43 true values
1241         $sdat4  = whereND $data, $mask4;   # $sdat4 is a [$n4,3,2] pdl
1242         $sdat43 = whereND $data, $mask43;  # $sdat43 is a [$n43,2] pdl
1243
1244       Just as with "where", you can use the returned value in an assignment.
1245       That means that both of these examples are valid:
1246
1247         # Used to create a new slice stored in $sdat4:
1248         $sdat4 = $data->whereND($mask4);
1249         $sdat4 .= 0;
1250         # Used in lvalue context:
1251         $data->whereND($mask4) .= 0;
1252
1253   whichND
1254       Return the coordinates of non-zero values in a mask.
1255
1256       WhichND returns the N-dimensional coordinates of each nonzero value in
1257       a mask PDL with any number of dimensions.  The returned values arrive
1258       as an array-of-vectors suitable for use in indexND or range.
1259
1260        $coords = whichND($mask);
1261
1262       returns a PDL containing the coordinates of the elements that are non-
1263       zero in $mask, suitable for use in indexND.  The 0th dimension contains
1264       the full coordinate listing of each point; the 1st dimension lists all
1265       the points.  For example, if $mask has rank 4 and 100 matching
1266       elements, then $coords has dimension 4x100.
1267
1268       If no such elements exist, then whichND returns a structured empty PDL:
1269       an Nx0 PDL that contains no values (but matches, threading-wise, with
1270       the vectors that would be produced if such elements existed).
1271
1272       DEPRECATED BEHAVIOR IN LIST CONTEXT:
1273
1274       whichND once delivered different values in list context than in scalar
1275       context, for historical reasons.  In list context, it returned the
1276       coordinates transposed, as a collection of 1-PDLs (one per dimension)
1277       in a list.  This usage is deprecated in PDL 2.4.10, and will cause a
1278       warning to be issued every time it is encountered.  To avoid the
1279       warning, you can set the global variable "$PDL::whichND" to 's' to get
1280       scalar behavior in all contexts, or to 'l' to get list behavior in list
1281       context.
1282
1283       In later versions of PDL, the deprecated behavior will disappear.
1284       Deprecated list context whichND expressions can be replaced with:
1285
1286           @list = $a->whichND->mv(0,-1)->dog;
1287
1288       SEE ALSO:
1289
1290       which finds coordinates of nonzero values in a 1-D mask.
1291
1292       where extracts values from a data PDL that are associated with nonzero
1293       values in a mask PDL.
1294
1295        pdl> $a=sequence(10,10,3,4)
1296        pdl> ($x, $y, $z, $w)=whichND($a == 203); p $x, $y, $z, $w
1297        [3] [0] [2] [0]
1298        pdl> print $a->at(list(cat($x,$y,$z,$w)))
1299        203
1300
1301   setops
1302       Implements simple set operations like union and intersection
1303
1304          Usage: $set = setops($a, <OPERATOR>, $b);
1305
1306       The operator can be "OR", "XOR" or "AND". This is then applied to $a
1307       viewed as a set and $b viewed as a set. Set theory says that a set may
1308       not have two or more identical elements, but setops takes care of this
1309       for you, so "$a=pdl(1,1,2)" is OK. The functioning is as follows:
1310
1311       "OR"
1312           The resulting vector will contain the elements that are either in
1313           $a or in $b or both. This is the union in set operation terms
1314
1315       "XOR"
1316           The resulting vector will contain the elements that are either in
1317           $a or $b, but not in both. This is
1318
1319                Union($a, $b) - Intersection($a, $b)
1320
1321           in set operation terms.
1322
1323       "AND"
1324           The resulting vector will contain the intersection of $a and $b, so
1325           the elements that are in both $a and $b. Note that for convenience
1326           this operation is also aliased to intersect.
1327
1328       It should be emphasized that these routines are used when one or both
1329       of the sets $a, $b are hard to calculate or that you get from a
1330       separate subroutine.
1331
1332       Finally IDL users might be familiar with Craig Markwardt's
1333       "cmset_op.pro" routine which has inspired this routine although it was
1334       written independently However the present routine has a few less
1335       options (but see the examples)
1336
1337       You will very often use these functions on an index vector, so that is
1338       what we will show here. We will in fact something slightly silly. First
1339       we will find all squares that are also cubes below 10000.
1340
1341       Create a sequence vector:
1342
1343         pdl> $x = sequence(10000)
1344
1345       Find all odd and even elements:
1346
1347         pdl> ($even, $odd) = which_both( ($x % 2) == 0)
1348
1349       Find all squares
1350
1351         pdl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
1352
1353       Find all cubes (being careful with roundoff error!)
1354
1355         pdl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
1356
1357       Then find all squares that are cubes:
1358
1359         pdl> $both = setops($squares, 'AND', $cubes)
1360
1361       And print these (assumes that "PDL::NiceSlice" is loaded!)
1362
1363         pdl> p $x($both)
1364          [0 1 64 729 4096]
1365
1366       Then find all numbers that are either cubes or squares, but not both:
1367
1368         pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)
1369
1370         pdl> p $cube_xor_square->nelem()
1371          112
1372
1373       So there are a total of 112 of these!
1374
1375       Finally find all odd squares:
1376
1377         pdl> $odd_squares = setops($squares, 'AND', $odd)
1378
1379       Another common occurrence is to want to get all objects that are in $a
1380       and in the complement of $b. But it is almost always best to create the
1381       complement explicitly since the universe that both are taken from is
1382       not known. Thus use which_both if possible to keep track of
1383       complements.
1384
1385       If this is impossible the best approach is to make a temporary:
1386
1387       This creates an index vector the size of the universe of the sets and
1388       set all elements in $b to 0
1389
1390         pdl> $tmp = ones($n_universe); $tmp($b) .= 0;
1391
1392       This then finds the complement of $b
1393
1394         pdl> $C_b = which($tmp == 1);
1395
1396       and this does the final selection:
1397
1398         pdl> $set = setops($a, 'AND', $C_b)
1399
1400   intersect
1401       Calculate the intersection of two piddles
1402
1403          Usage: $set = intersect($a, $b);
1404
1405       This routine is merely a simple interface to setops. See that for more
1406       information
1407
1408       Find all numbers less that 100 that are of the form 2*y and 3*x
1409
1410        pdl> $x=sequence(100)
1411        pdl> $factor2 = which( ($x % 2) == 0)
1412        pdl> $factor3 = which( ($x % 3) == 0)
1413        pdl> $ii=intersect($factor2, $factor3)
1414        pdl> p $x($ii)
1415        [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
1416

AUTHOR

1418       Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu).
1419       Contributions by Christian Soeller (c.soeller@auckland.ac.nz), Karl
1420       Glazebrook (kgb@aaoepp.aao.gov.au), Craig DeForest
1421       (deforest@boulder.swri.edu) and Jarle Brinchmann (jarle@astro.up.pt)
1422       All rights reserved. There is no warranty. You are allowed to
1423       redistribute this software / documentation under certain conditions.
1424       For details, see the file COPYING in the PDL distribution. If this file
1425       is separated from the PDL distribution, the copyright notice should be
1426       included in the file.
1427
1428       Updated for CPAN viewing compatibility by David Mertens.
1429
1430
1431
1432perl v5.30.0                      2019-09-05                      Primitive(3)
Impressum