"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/GENERATED/PDL/Primitive.pm" (28 May 2022, 94414 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Primitive.pm" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.079_vs_2.080.

A hint: This file contains one or more very long lines, so maybe it is better readable using the pure text view mode that shows the contents as wrapped lines within the browser window.


    1 #
    2 # GENERATED WITH PDL::PP! Don't modify!
    3 #
    4 package PDL::Primitive;
    5 
    6 our @EXPORT_OK = qw(inner outer matmult innerwt inner2 inner2d inner2t crossp norm indadd conv1d in uniq uniqind uniqvec hclip lclip clip clip wtstat statsover stats histogram whistogram histogram2d whistogram2d fibonacci append axisvalues cmpvec eqvec enumvec enumvecg vsearchvec unionvec intersectvec setdiffvec union_sorted intersect_sorted setdiff_sorted srand random randsym grandom vsearch vsearch_sample vsearch_insert_leftmost vsearch_insert_rightmost vsearch_match vsearch_bin_inclusive vsearch_bin_exclusive interpolate interpol interpND one2nd which which_both where where_both whereND whichND setops intersect );
    7 our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
    8 
    9 use PDL::Core;
   10 use PDL::Exporter;
   11 use DynaLoader;
   12 
   13 
   14    
   15    our @ISA = ( 'PDL::Exporter','DynaLoader' );
   16    push @PDL::Core::PP, __PACKAGE__;
   17    bootstrap PDL::Primitive ;
   18 
   19 
   20 
   21 
   22 
   23 
   24 #line 6 "primitive.pd"
   25 
   26 use strict;
   27 use warnings;
   28 use PDL::Slices;
   29 use Carp;
   30 
   31 { package PDL;
   32   use overload (
   33     'x' => sub {
   34       PDL::Primitive::matmult(@_[0,1], my $foo=$_[0]->null());
   35       $foo;
   36     },
   37   );
   38 }
   39 
   40 =head1 NAME
   41 
   42 PDL::Primitive - primitive operations for pdl
   43 
   44 =head1 DESCRIPTION
   45 
   46 This module provides some primitive and useful functions defined
   47 using PDL::PP and able to use the new indexing tricks.
   48 
   49 See L<PDL::Indexing> for how to use indices creatively.
   50 For explanation of the signature format, see L<PDL::PP>.
   51 
   52 =head1 SYNOPSIS
   53 
   54  # Pulls in PDL::Primitive, among other modules.
   55  use PDL;
   56 
   57  # Only pull in PDL::Primitive:
   58  use PDL::Primitive;
   59 
   60 =cut
   61 #line 62 "Primitive.pm"
   62 
   63 
   64 
   65 
   66 
   67 
   68 =head1 FUNCTIONS
   69 
   70 =cut
   71 
   72 
   73 
   74 
   75 #line 948 "../../blib/lib/PDL/PP.pm"
   76 
   77 
   78 
   79 =head2 inner
   80 
   81 =for sig
   82 
   83   Signature: (a(n); b(n); [o]c())
   84 
   85 
   86 =for ref
   87 
   88 Inner product over one dimension
   89 
   90  c = sum_i a_i * b_i
   91 
   92 
   93 =for bad
   94 
   95 =for bad
   96 
   97 If C<a() * b()> contains only bad data,
   98 C<c()> is set bad. Otherwise C<c()> will have its bad flag cleared,
   99 as it will not contain any bad values.
  100 
  101 
  102 =cut
  103 #line 104 "Primitive.pm"
  104 
  105 
  106 
  107 #line 950 "../../blib/lib/PDL/PP.pm"
  108 
  109 *inner = \&PDL::inner;
  110 #line 111 "Primitive.pm"
  111 
  112 
  113 
  114 #line 948 "../../blib/lib/PDL/PP.pm"
  115 
  116 
  117 
  118 =head2 outer
  119 
  120 =for sig
  121 
  122   Signature: (a(n); b(m); [o]c(n,m))
  123 
  124 
  125 =for ref
  126 
  127 outer product over one dimension
  128 
  129 Naturally, it is possible to achieve the effects of outer
  130 product simply by broadcasting over the "C<*>"
  131 operator but this function is provided for convenience.
  132 
  133 
  134 =for bad
  135 
  136 outer processes bad values.
  137 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  138 
  139 
  140 =cut
  141 #line 142 "Primitive.pm"
  142 
  143 
  144 
  145 #line 950 "../../blib/lib/PDL/PP.pm"
  146 
  147 *outer = \&PDL::outer;
  148 #line 149 "Primitive.pm"
  149 
  150 
  151 
  152 #line 104 "primitive.pd"
  153 
  154 =head2 x
  155 
  156 =for sig
  157 
  158  Signature: (a(i,z), b(x,i),[o]c(x,z))
  159 
  160 =for ref
  161 
  162 Matrix multiplication
  163 
  164 PDL overloads the C<x> operator (normally the repeat operator) for
  165 matrix multiplication.  The number of columns (size of the 0
  166 dimension) in the left-hand argument must normally equal the number of
  167 rows (size of the 1 dimension) in the right-hand argument.
  168 
  169 Row vectors are represented as (N x 1) two-dimensional PDLs, or you
  170 may be sloppy and use a one-dimensional PDL.  Column vectors are
  171 represented as (1 x N) two-dimensional PDLs.
  172 
  173 Broadcasting occurs in the usual way, but as both the 0 and 1 dimension
  174 (if present) are included in the operation, you must be sure that
  175 you don't try to broadcast over either of those dims.
  176 
  177 Of note, due to how Perl v5.14.0 and above implement operator overloading of
  178 the C<x> operator, the use of parentheses for the left operand creates a list
  179 context, that is
  180 
  181  pdl> ( $x * $y ) x $z
  182  ERROR: Argument "..." isn't numeric in repeat (x) ...
  183 
  184 treats C<$z> as a numeric count for the list repeat operation and does not call
  185 the scalar form of the overloaded operator. To use the operator in this case,
  186 use a scalar context:
  187 
  188  pdl> scalar( $x * $y ) x $z
  189 
  190 or by calling L</matmult> directly:
  191 
  192  pdl> ( $x * $y )->matmult( $z )
  193 
  194 EXAMPLES
  195 
  196 Here are some simple ways to define vectors and matrices:
  197 
  198  pdl> $r = pdl(1,2);                # A row vector
  199  pdl> $c = pdl([[3],[4]]);          # A column vector
  200  pdl> $c = pdl(3,4)->(*1);          # A column vector, using NiceSlice
  201  pdl> $m = pdl([[1,2],[3,4]]);      # A 2x2 matrix
  202 
  203 Now that we have a few objects prepared, here is how to
  204 matrix-multiply them:
  205 
  206  pdl> print $r x $m                 # row x matrix = row
  207  [
  208   [ 7 10]
  209  ]
  210 
  211  pdl> print $m x $r                 # matrix x row = ERROR
  212  PDL: Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1
  213 
  214  pdl> print $m x $c                 # matrix x column = column
  215  [
  216   [ 5]
  217   [11]
  218  ]
  219 
  220  pdl> print $m x 2                  # Trivial case: scalar mult.
  221  [
  222   [2 4]
  223   [6 8]
  224  ]
  225 
  226  pdl> print $r x $c                 # row x column = scalar
  227  [
  228   [11]
  229  ]
  230 
  231  pdl> print $c x $r                 # column x row = matrix
  232  [
  233   [3 6]
  234   [4 8]
  235  ]
  236 
  237 INTERNALS
  238 
  239 The mechanics of the multiplication are carried out by the
  240 L</matmult> method.
  241 
  242 =cut
  243 #line 244 "Primitive.pm"
  244 
  245 
  246 
  247 #line 948 "../../blib/lib/PDL/PP.pm"
  248 
  249 
  250 
  251 =head2 matmult
  252 
  253 =for sig
  254 
  255   Signature: (a(t,h); b(w,t); [o]c(w,h))
  256 
  257 =for ref
  258 
  259 Matrix multiplication
  260 
  261 Notionally, matrix multiplication $x x $y is equivalent to the
  262 broadcasting expression
  263 
  264     $x->dummy(1)->inner($y->xchg(0,1)->dummy(2),$c);
  265 
  266 but for large matrices that breaks CPU cache and is slow.  Instead,
  267 matmult calculates its result in 32x32x32 tiles, to keep the memory
  268 footprint within cache as long as possible on most modern CPUs.
  269 
  270 For usage, see L</x>, a description of the overloaded 'x' operator
  271 
  272 
  273 
  274 =for bad
  275 
  276 matmult ignores the bad-value flag of the input ndarrays.
  277 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  278 
  279 
  280 =cut
  281 #line 282 "Primitive.pm"
  282 
  283 
  284 
  285 #line 949 "../../blib/lib/PDL/PP.pm"
  286 
  287 sub PDL::matmult {
  288     my ($x,$y,$c) = @_;
  289     $y = PDL->topdl($y);
  290     $c = PDL->null unless do { local $@; eval { $c->isa('PDL') } };
  291     while($x->getndims < 2) {$x = $x->dummy(-1)}
  292     while($y->getndims < 2) {$y = $y->dummy(-1)}
  293     return ($c .= $x * $y) if( ($x->dim(0)==1 && $x->dim(1)==1) ||
  294                                ($y->dim(0)==1 && $y->dim(1)==1) );
  295     barf sprintf 'Dim mismatch in matmult of [%1$dx%2$d] x [%3$dx%4$d]: %1$d != %4$d',$x->dim(0),$x->dim(1),$y->dim(0),$y->dim(1)
  296       if $y->dim(1) != $x->dim(0);
  297     PDL::_matmult_int($x,$y,$c);
  298     $c;
  299 }
  300 #line 301 "Primitive.pm"
  301 
  302 
  303 
  304 #line 950 "../../blib/lib/PDL/PP.pm"
  305 
  306 *matmult = \&PDL::matmult;
  307 #line 308 "Primitive.pm"
  308 
  309 
  310 
  311 #line 948 "../../blib/lib/PDL/PP.pm"
  312 
  313 
  314 
  315 =head2 innerwt
  316 
  317 =for sig
  318 
  319   Signature: (a(n); b(n); c(n); [o]d())
  320 
  321 
  322 
  323 =for ref
  324 
  325 Weighted (i.e. triple) inner product
  326 
  327  d = sum_i a(i) b(i) c(i)
  328 
  329 
  330 
  331 =for bad
  332 
  333 innerwt processes bad values.
  334 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  335 
  336 
  337 =cut
  338 #line 339 "Primitive.pm"
  339 
  340 
  341 
  342 #line 950 "../../blib/lib/PDL/PP.pm"
  343 
  344 *innerwt = \&PDL::innerwt;
  345 #line 346 "Primitive.pm"
  346 
  347 
  348 
  349 #line 948 "../../blib/lib/PDL/PP.pm"
  350 
  351 
  352 
  353 =head2 inner2
  354 
  355 =for sig
  356 
  357   Signature: (a(n); b(n,m); c(m); [o]d())
  358 
  359 
  360 =for ref
  361 
  362 Inner product of two vectors and a matrix
  363 
  364  d = sum_ij a(i) b(i,j) c(j)
  365 
  366 Note that you should probably not broadcast over C<a> and C<c> since that would be
  367 very wasteful. Instead, you should use a temporary for C<b*c>.
  368 
  369 
  370 =for bad
  371 
  372 inner2 processes bad values.
  373 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  374 
  375 
  376 =cut
  377 #line 378 "Primitive.pm"
  378 
  379 
  380 
  381 #line 950 "../../blib/lib/PDL/PP.pm"
  382 
  383 *inner2 = \&PDL::inner2;
  384 #line 385 "Primitive.pm"
  385 
  386 
  387 
  388 #line 948 "../../blib/lib/PDL/PP.pm"
  389 
  390 
  391 
  392 =head2 inner2d
  393 
  394 =for sig
  395 
  396   Signature: (a(n,m); b(n,m); [o]c())
  397 
  398 
  399 
  400 =for ref
  401 
  402 Inner product over 2 dimensions.
  403 
  404 Equivalent to
  405 
  406  $c = inner($x->clump(2), $y->clump(2))
  407 
  408 
  409 
  410 =for bad
  411 
  412 inner2d processes bad values.
  413 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  414 
  415 
  416 =cut
  417 #line 418 "Primitive.pm"
  418 
  419 
  420 
  421 #line 950 "../../blib/lib/PDL/PP.pm"
  422 
  423 *inner2d = \&PDL::inner2d;
  424 #line 425 "Primitive.pm"
  425 
  426 
  427 
  428 #line 948 "../../blib/lib/PDL/PP.pm"
  429 
  430 
  431 
  432 =head2 inner2t
  433 
  434 =for sig
  435 
  436   Signature: (a(j,n); b(n,m); c(m,k); [t]tmp(n,k); [o]d(j,k)))
  437 
  438 
  439 =for ref
  440 
  441 Efficient Triple matrix product C<a*b*c>
  442 
  443 Efficiency comes from by using the temporary C<tmp>. This operation only
  444 scales as C<N**3> whereas broadcasting using L</inner2> would scale
  445 as C<N**4>.
  446 
  447 The reason for having this routine is that you do not need to
  448 have the same broadcast-dimensions for C<tmp> as for the other arguments,
  449 which in case of large numbers of matrices makes this much more
  450 memory-efficient.
  451 
  452 It is hoped that things like this could be taken care of as a kind of
  453 closures at some point.
  454 
  455 
  456 
  457 =for bad
  458 
  459 inner2t processes bad values.
  460 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  461 
  462 
  463 =cut
  464 #line 465 "Primitive.pm"
  465 
  466 
  467 
  468 #line 950 "../../blib/lib/PDL/PP.pm"
  469 
  470 *inner2t = \&PDL::inner2t;
  471 #line 472 "Primitive.pm"
  472 
  473 
  474 
  475 #line 948 "../../blib/lib/PDL/PP.pm"
  476 
  477 
  478 
  479 =head2 crossp
  480 
  481 =for sig
  482 
  483   Signature: (a(tri=3); b(tri); [o] c(tri))
  484 
  485 
  486 =for ref
  487 
  488 Cross product of two 3D vectors
  489 
  490 After
  491 
  492 =for example
  493 
  494  $c = crossp $x, $y
  495 
  496 the inner product C<$c*$x> and C<$c*$y> will be zero, i.e. C<$c> is
  497 orthogonal to C<$x> and C<$y>
  498 
  499 
  500 
  501 =for bad
  502 
  503 crossp does not process bad values.
  504 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  505 
  506 
  507 =cut
  508 #line 509 "Primitive.pm"
  509 
  510 
  511 
  512 #line 950 "../../blib/lib/PDL/PP.pm"
  513 
  514 *crossp = \&PDL::crossp;
  515 #line 516 "Primitive.pm"
  516 
  517 
  518 
  519 #line 948 "../../blib/lib/PDL/PP.pm"
  520 
  521 
  522 
  523 =head2 norm
  524 
  525 =for sig
  526 
  527   Signature: (vec(n); [o] norm(n))
  528 
  529 =for ref
  530 
  531 Normalises a vector to unit Euclidean length
  532 
  533 =for bad
  534 
  535 norm processes bad values.
  536 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  537 
  538 
  539 =cut
  540 #line 541 "Primitive.pm"
  541 
  542 
  543 
  544 #line 950 "../../blib/lib/PDL/PP.pm"
  545 
  546 *norm = \&PDL::norm;
  547 #line 548 "Primitive.pm"
  548 
  549 
  550 
  551 #line 948 "../../blib/lib/PDL/PP.pm"
  552 
  553 
  554 
  555 =head2 indadd
  556 
  557 =for sig
  558 
  559   Signature: (input(n); indx ind(n); [io] sum(m))
  560 
  561 
  562 =for ref
  563 
  564 Broadcasting index add: add C<input> to the C<ind> element of C<sum>, i.e:
  565 
  566  sum(ind) += input
  567 
  568 =for example
  569 
  570 Simple example:
  571 
  572   $x = 2;
  573   $ind = 3;
  574   $sum = zeroes(10);
  575   indadd($x,$ind, $sum);
  576   print $sum
  577   #Result: ( 2 added to element 3 of $sum)
  578   # [0 0 0 2 0 0 0 0 0 0]
  579 
  580 Broadcasting example:
  581 
  582   $x = pdl( 1,2,3);
  583   $ind = pdl( 1,4,6);
  584   $sum = zeroes(10);
  585   indadd($x,$ind, $sum);
  586   print $sum."\n";
  587   #Result: ( 1, 2, and 3 added to elements 1,4,6 $sum)
  588   # [0 1 0 0 2 0 3 0 0 0]
  589 
  590 
  591 
  592 =for bad
  593 
  594 =for bad
  595 
  596 The routine barfs on bad indices, and bad inputs set target outputs bad.
  597 
  598 
  599 
  600 =cut
  601 #line 602 "Primitive.pm"
  602 
  603 
  604 
  605 #line 950 "../../blib/lib/PDL/PP.pm"
  606 
  607 *indadd = \&PDL::indadd;
  608 #line 609 "Primitive.pm"
  609 
  610 
  611 
  612 #line 948 "../../blib/lib/PDL/PP.pm"
  613 
  614 
  615 
  616 =head2 conv1d
  617 
  618 =for sig
  619 
  620   Signature: (a(m); kern(p); [o]b(m); int reflect)
  621 
  622 
  623 =for ref
  624 
  625 1D convolution along first dimension
  626 
  627 The m-th element of the discrete convolution of an input ndarray
  628 C<$a> of size C<$M>, and a kernel ndarray C<$kern> of size C<$P>, is
  629 calculated as
  630 
  631                               n = ($P-1)/2
  632                               ====
  633                               \
  634   ($a conv1d $kern)[m]   =     >      $a_ext[m - n] * $kern[n]
  635                               /
  636                               ====
  637                               n = -($P-1)/2
  638 
  639 where C<$a_ext> is either the periodic (or reflected) extension of
  640 C<$a> so it is equal to C<$a> on C< 0..$M-1 > and equal to the
  641 corresponding periodic/reflected image of C<$a> outside that range.
  642 
  643 
  644 =for example
  645 
  646   $con = conv1d sequence(10), pdl(-1,0,1);
  647 
  648   $con = conv1d sequence(10), pdl(-1,0,1), {Boundary => 'reflect'};
  649 
  650 By default, periodic boundary conditions are assumed (i.e. wrap around).
  651 Alternatively, you can request reflective boundary conditions using
  652 the C<Boundary> option:
  653 
  654   {Boundary => 'reflect'} # case in 'reflect' doesn't matter
  655 
  656 The convolution is performed along the first dimension. To apply it across
  657 another dimension use the slicing routines, e.g.
  658 
  659   $y = $x->mv(2,0)->conv1d($kernel)->mv(0,2); # along third dim
  660 
  661 This function is useful for broadcasted filtering of 1D signals.
  662 
  663 Compare also L<conv2d|PDL::Image2D/conv2d>, L<convolve|PDL::ImageND/convolve>,
  664 L<fftconvolve|PDL::FFT/fftconvolve()>, L<fftwconv|PDL::FFTW/fftwconv>,
  665 L<rfftwconv|PDL::FFTW/rfftwconv>
  666 
  667 =for bad
  668 
  669 WARNING: C<conv1d> processes bad values in its inputs as
  670 the numeric value of C<< $pdl->badvalue >> so it is not
  671 recommended for processing pdls with bad values in them
  672 unless special care is taken.
  673 
  674 
  675 
  676 =for bad
  677 
  678 conv1d ignores the bad-value flag of the input ndarrays.
  679 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  680 
  681 
  682 =cut
  683 #line 684 "Primitive.pm"
  684 
  685 
  686 
  687 #line 949 "../../blib/lib/PDL/PP.pm"
  688 
  689 
  690 
  691 sub PDL::conv1d {
  692    my $opt = pop @_ if ref($_[-1]) eq 'HASH';
  693    die 'Usage: conv1d( a(m), kern(p), [o]b(m), {Options} )'
  694       if @_<2 || @_>3;
  695    my($x,$kern) = @_;
  696    my $c = @_ == 3 ? $_[2] : PDL->null;
  697    PDL::_conv1d_int($x,$kern,$c,
  698              !(defined $opt && exists $$opt{Boundary}) ? 0 :
  699              lc $$opt{Boundary} eq "reflect");
  700    return $c;
  701 }
  702 #line 703 "Primitive.pm"
  703 
  704 
  705 
  706 #line 950 "../../blib/lib/PDL/PP.pm"
  707 
  708 *conv1d = \&PDL::conv1d;
  709 #line 710 "Primitive.pm"
  710 
  711 
  712 
  713 #line 948 "../../blib/lib/PDL/PP.pm"
  714 
  715 
  716 
  717 =head2 in
  718 
  719 =for sig
  720 
  721   Signature: (a(); b(n); [o] c())
  722 
  723 
  724 =for ref
  725 
  726 test if a is in the set of values b
  727 
  728 =for example
  729 
  730    $goodmsk = $labels->in($goodlabels);
  731    print pdl(3,1,4,6,2)->in(pdl(2,3,3));
  732   [1 0 0 0 1]
  733 
  734 C<in> is akin to the I<is an element of> of set theory. In principle,
  735 PDL broadcasting could be used to achieve its functionality by using a
  736 construct like
  737 
  738    $msk = ($labels->dummy(0) == $goodlabels)->orover;
  739 
  740 However, C<in> doesn't create a (potentially large) intermediate
  741 and is generally faster.
  742 
  743 
  744 
  745 =for bad
  746 
  747 in does not process bad values.
  748 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  749 
  750 
  751 =cut
  752 #line 753 "Primitive.pm"
  753 
  754 
  755 
  756 #line 950 "../../blib/lib/PDL/PP.pm"
  757 
  758 *in = \&PDL::in;
  759 #line 760 "Primitive.pm"
  760 
  761 
  762 
  763 #line 697 "primitive.pd"
  764 
  765 =head2 uniq
  766 
  767 =for ref
  768 
  769 return all unique elements of an ndarray
  770 
  771 The unique elements are returned in ascending order.
  772 
  773 =for example
  774 
  775   PDL> p pdl(2,2,2,4,0,-1,6,6)->uniq
  776   [-1 0 2 4 6]     # 0 is returned 2nd (sorted order)
  777 
  778   PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniq
  779   [-1 2 4 6 nan]   # NaN value is returned at end
  780 
  781 Note: The returned pdl is 1D; any structure of the input
  782 ndarray is lost.  C<NaN> values are never compare equal to
  783 any other values, even themselves.  As a result, they are
  784 always unique. C<uniq> returns the NaN values at the end
  785 of the result ndarray.  This follows the Matlab usage.
  786 
  787 See L</uniqind> if you need the indices of the unique
  788 elements rather than the values.
  789 
  790 =for bad
  791 
  792 Bad values are not considered unique by uniq and are ignored.
  793 
  794  $x=sequence(10);
  795  $x=$x->setbadif($x%3);
  796  print $x->uniq;
  797  [0 3 6 9]
  798 
  799 =cut
  800 
  801 *uniq = \&PDL::uniq;
  802 # return unique elements of array
  803 # find as jumps in the sorted array
  804 # flattens in the process
  805 sub PDL::uniq {
  806    use PDL::Core 'barf';
  807    my ($arr) = @_;
  808    return $arr if($arr->nelem == 0); # The null list is unique (CED)
  809    my $srt  = $arr->clump(-1)->where($arr==$arr)->qsort;  # no NaNs or BADs for qsort
  810    my $nans = $arr->clump(-1)->where($arr!=$arr);
  811    my $uniq = ($srt->nelem > 0) ? $srt->where($srt != $srt->rotate(-1)) : $srt;
  812    # make sure we return something if there is only one value
  813    my $answ = $nans;  # NaN values always uniq
  814    if ( $uniq->nelem > 0 ) {
  815       $answ = $uniq->append($answ);
  816    } else {
  817       $answ = ( ($srt->nelem == 0) ?  $srt : PDL::pdl( ref($srt), [$srt->index(0)] ) )->append($answ);
  818    }
  819    return $answ;
  820 }
  821 #line 822 "Primitive.pm"
  822 
  823 
  824 
  825 #line 757 "primitive.pd"
  826 
  827 =head2 uniqind
  828 
  829 =for ref
  830 
  831 Return the indices of all unique elements of an ndarray
  832 The order is in the order of the values to be consistent
  833 with uniq. C<NaN> values never compare equal with any
  834 other value and so are always unique.  This follows the
  835 Matlab usage.
  836 
  837 =for example
  838 
  839   PDL> p pdl(2,2,2,4,0,-1,6,6)->uniqind
  840   [5 4 1 3 6]     # the 0 at index 4 is returned 2nd, but...
  841 
  842   PDL> p pdl(2,2,2,4,nan,-1,6,6)->uniqind
  843   [5 1 3 6 4]     # ...the NaN at index 4 is returned at end
  844 
  845 
  846 Note: The returned pdl is 1D; any structure of the input
  847 ndarray is lost.
  848 
  849 See L</uniq> if you want the unique values instead of the
  850 indices.
  851 
  852 =for bad
  853 
  854 Bad values are not considered unique by uniqind and are ignored.
  855 
  856 =cut
  857 
  858 *uniqind = \&PDL::uniqind;
  859 # return unique elements of array
  860 # find as jumps in the sorted array
  861 # flattens in the process
  862 sub PDL::uniqind {
  863   use PDL::Core 'barf';
  864   my ($arr) = @_;
  865   return $arr if($arr->nelem == 0); # The null list is unique (CED)
  866   # Different from uniq we sort and store the result in an intermediary
  867   my $aflat = $arr->flat;
  868   my $nanind = which($aflat!=$aflat);                        # NaN indexes
  869   my $good = PDL->sequence(indx, $aflat->dims)->where($aflat==$aflat);  # good indexes
  870   my $i_srt = $aflat->where($aflat==$aflat)->qsorti;         # no BAD or NaN values for qsorti
  871   my $srt = $aflat->where($aflat==$aflat)->index($i_srt);
  872   my $uniqind;
  873   if ($srt->nelem > 0) {
  874      $uniqind = which($srt != $srt->rotate(-1));
  875      $uniqind = $i_srt->slice('0') if $uniqind->isempty;
  876   } else {
  877      $uniqind = which($srt);
  878   }
  879   # Now map back to the original space
  880   my $ansind = $nanind;
  881   if ( $uniqind->nelem > 0 ) {
  882      $ansind = ($good->index($i_srt->index($uniqind)))->append($ansind);
  883   } else {
  884      $ansind = $uniqind->append($ansind);
  885   }
  886   return $ansind;
  887 }
  888 #line 889 "Primitive.pm"
  889 
  890 
  891 
  892 #line 823 "primitive.pd"
  893 
  894 =head2 uniqvec
  895 
  896 =for ref
  897 
  898 Return all unique vectors out of a collection
  899 
  900   NOTE: If any vectors in the input ndarray have NaN values
  901   they are returned at the end of the non-NaN ones.  This is
  902   because, by definition, NaN values never compare equal with
  903   any other value.
  904 
  905   NOTE: The current implementation does not sort the vectors
  906   containing NaN values.
  907 
  908 The unique vectors are returned in lexicographically sorted
  909 ascending order. The 0th dimension of the input PDL is treated
  910 as a dimensional index within each vector, and the 1st and any
  911 higher dimensions are taken to run across vectors. The return
  912 value is always 2D; any structure of the input PDL (beyond using
  913 the 0th dimension for vector index) is lost.
  914 
  915 See also L</uniq> for a unique list of scalars; and
  916 L<qsortvec|PDL::Ufunc/qsortvec> for sorting a list of vectors
  917 lexicographcally.
  918 
  919 =for bad
  920 
  921 If a vector contains all bad values, it is ignored as in L</uniq>.
  922 If some of the values are good, it is treated as a normal vector. For
  923 example, [1 2 BAD] and [BAD 2 3] could be returned, but [BAD BAD BAD]
  924 could not.  Vectors containing BAD values will be returned after any
  925 non-NaN and non-BAD containing vectors, followed by the NaN vectors.
  926 
  927 =cut
  928 
  929 sub PDL::uniqvec {
  930    my($pdl) = shift;
  931 
  932    return $pdl if ( $pdl->nelem == 0 || $pdl->ndims < 2 );
  933    return $pdl if ( $pdl->slice("(0)")->nelem < 2 );                     # slice isn't cheap but uniqvec isn't either
  934 
  935    my $pdl2d = $pdl->clump(1..$pdl->ndims-1);
  936    my $ngood = $pdl2d->ngoodover;
  937    $pdl2d = $pdl2d->mv(0,-1)->dice($ngood->which)->mv(-1,0);             # remove all-BAD vectors
  938 
  939    my $numnan = ($pdl2d!=$pdl2d)->sumover;                                  # works since no all-BADs to confuse
  940 
  941    my $presrt = $pdl2d->mv(0,-1)->dice($numnan->not->which)->mv(0,-1);      # remove vectors with any NaN values
  942    my $nanvec = $pdl2d->mv(0,-1)->dice($numnan->which)->mv(0,-1);           # the vectors with any NaN values
  943 
  944    my $srt = $presrt->qsortvec->mv(0,-1);                                   # BADs are sorted by qsortvec
  945    my $srtdice = $srt;
  946    my $somebad = null;
  947    if ($srt->badflag) {
  948       $srtdice = $srt->dice($srt->mv(0,-1)->nbadover->not->which);
  949       $somebad = $srt->dice($srt->mv(0,-1)->nbadover->which);
  950    }
  951 
  952    my $uniq = $srtdice->nelem > 0
  953      ? ($srtdice != $srtdice->rotate(-1))->mv(0,-1)->orover->which
  954      : $srtdice->orover->which;
  955 
  956    my $ans = $uniq->nelem > 0 ? $srtdice->dice($uniq) :
  957       ($srtdice->nelem > 0) ? $srtdice->slice("0,:") :
  958       $srtdice;
  959    return $ans->append($somebad)->append($nanvec->mv(0,-1))->mv(0,-1);
  960 }
  961 #line 962 "Primitive.pm"
  962 
  963 
  964 
  965 #line 948 "../../blib/lib/PDL/PP.pm"
  966 
  967 
  968 
  969 =head2 hclip
  970 
  971 =for sig
  972 
  973   Signature: (a(); b(); [o] c())
  974 
  975 =for ref
  976 
  977 clip (threshold) C<$a> by C<$b> (C<$b> is upper bound)
  978 
  979 =for bad
  980 
  981 hclip processes bad values.
  982 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  983 
  984 
  985 =cut
  986 #line 987 "Primitive.pm"
  987 
  988 
  989 
  990 #line 949 "../../blib/lib/PDL/PP.pm"
  991 
  992 sub PDL::hclip {
  993    my ($x,$y) = @_;
  994    my $c;
  995    if ($x->is_inplace) {
  996        $x->set_inplace(0); $c = $x;
  997    } elsif (@_ > 2) {$c=$_[2]} else {$c=PDL->nullcreate($x)}
  998    PDL::_hclip_int($x,$y,$c);
  999    return $c;
 1000 }
 1001 #line 1002 "Primitive.pm"
 1002 
 1003 
 1004 
 1005 #line 950 "../../blib/lib/PDL/PP.pm"
 1006 
 1007 *hclip = \&PDL::hclip;
 1008 #line 1009 "Primitive.pm"
 1009 
 1010 
 1011 
 1012 #line 948 "../../blib/lib/PDL/PP.pm"
 1013 
 1014 
 1015 
 1016 =head2 lclip
 1017 
 1018 =for sig
 1019 
 1020   Signature: (a(); b(); [o] c())
 1021 
 1022 =for ref
 1023 
 1024 clip (threshold) C<$a> by C<$b> (C<$b> is lower bound)
 1025 
 1026 =for bad
 1027 
 1028 lclip processes bad values.
 1029 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1030 
 1031 
 1032 =cut
 1033 #line 1034 "Primitive.pm"
 1034 
 1035 
 1036 
 1037 #line 949 "../../blib/lib/PDL/PP.pm"
 1038 
 1039 sub PDL::lclip {
 1040    my ($x,$y) = @_;
 1041    my $c;
 1042    if ($x->is_inplace) {
 1043        $x->set_inplace(0); $c = $x;
 1044    } elsif (@_ > 2) {$c=$_[2]} else {$c=PDL->nullcreate($x)}
 1045    PDL::_lclip_int($x,$y,$c);
 1046    return $c;
 1047 }
 1048 #line 1049 "Primitive.pm"
 1049 
 1050 
 1051 
 1052 #line 950 "../../blib/lib/PDL/PP.pm"
 1053 
 1054 *lclip = \&PDL::lclip;
 1055 #line 1056 "Primitive.pm"
 1056 
 1057 
 1058 
 1059 #line 933 "primitive.pd"
 1060 
 1061 =head2 clip
 1062 
 1063 =for ref
 1064 
 1065 Clip (threshold) an ndarray by (optional) upper or lower bounds.
 1066 
 1067 =for usage
 1068 
 1069  $y = $x->clip(0,3);
 1070  $c = $x->clip(undef, $x);
 1071 
 1072 =for bad
 1073 
 1074 clip handles bad values since it is just a
 1075 wrapper around L</hclip> and
 1076 L</lclip>.
 1077 
 1078 =cut
 1079 #line 1080 "Primitive.pm"
 1080 
 1081 
 1082 
 1083 #line 948 "../../blib/lib/PDL/PP.pm"
 1084 
 1085 
 1086 
 1087 =head2 clip
 1088 
 1089 =for sig
 1090 
 1091   Signature: (a(); l(); h(); [o] c())
 1092 
 1093 
 1094 =for ref
 1095 
 1096 info not available
 1097 
 1098 
 1099 =for bad
 1100 
 1101 clip processes bad values.
 1102 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1103 
 1104 
 1105 =cut
 1106 #line 1107 "Primitive.pm"
 1107 
 1108 
 1109 
 1110 #line 949 "../../blib/lib/PDL/PP.pm"
 1111 
 1112 *clip = \&PDL::clip;
 1113 sub PDL::clip {
 1114   my($x, $l, $h) = @_;
 1115   my $d;
 1116   unless(defined($l) || defined($h)) {
 1117       # Deal with pathological case
 1118       if($x->is_inplace) {
 1119       $x->set_inplace(0);
 1120       return $x;
 1121       } else {
 1122       return $x->copy;
 1123       }
 1124   }
 1125 
 1126   if($x->is_inplace) {
 1127       $x->set_inplace(0); $d = $x
 1128   } elsif (@_ > 3) {
 1129       $d=$_[3]
 1130   } else {
 1131       $d = PDL->nullcreate($x);
 1132   }
 1133   if(defined($l) && defined($h)) {
 1134       PDL::_clip_int($x,$l,$h,$d);
 1135   } elsif( defined($l) ) {
 1136       PDL::_lclip_int($x,$l,$d);
 1137   } elsif( defined($h) ) {
 1138       PDL::_hclip_int($x,$h,$d);
 1139   } else {
 1140       die "This can't happen (clip contingency) - file a bug";
 1141   }
 1142 
 1143   return $d;
 1144 }
 1145 #line 1146 "Primitive.pm"
 1146 
 1147 
 1148 
 1149 #line 950 "../../blib/lib/PDL/PP.pm"
 1150 
 1151 *clip = \&PDL::clip;
 1152 #line 1153 "Primitive.pm"
 1153 
 1154 
 1155 
 1156 #line 948 "../../blib/lib/PDL/PP.pm"
 1157 
 1158 
 1159 
 1160 =head2 wtstat
 1161 
 1162 =for sig
 1163 
 1164   Signature: (a(n); wt(n); avg(); [o]b(); int deg)
 1165 
 1166 
 1167 =for ref
 1168 
 1169 Weighted statistical moment of given degree
 1170 
 1171 This calculates a weighted statistic over the vector C<a>.
 1172 The formula is
 1173 
 1174  b() = (sum_i wt_i * (a_i ** degree - avg)) / (sum_i wt_i)
 1175 
 1176 
 1177 =for bad
 1178 
 1179 =for bad
 1180 
 1181 Bad values are ignored in any calculation; C<$b> will only
 1182 have its bad flag set if the output contains any bad data.
 1183 
 1184 
 1185 =cut
 1186 #line 1187 "Primitive.pm"
 1187 
 1188 
 1189 
 1190 #line 950 "../../blib/lib/PDL/PP.pm"
 1191 
 1192 *wtstat = \&PDL::wtstat;
 1193 #line 1194 "Primitive.pm"
 1194 
 1195 
 1196 
 1197 #line 948 "../../blib/lib/PDL/PP.pm"
 1198 
 1199 
 1200 
 1201 =head2 statsover
 1202 
 1203 =for sig
 1204 
 1205   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())
 1206 
 1207 
 1208 =for ref
 1209 
 1210 Calculate useful statistics over a dimension of an ndarray
 1211 
 1212 =for usage
 1213 
 1214   ($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($ndarray, $weights);
 1215 
 1216 This utility function calculates various useful
 1217 quantities of an ndarray. These are:
 1218 
 1219 =over 3
 1220 
 1221 =item * the mean:
 1222 
 1223   MEAN = sum (x)/ N
 1224 
 1225 with C<N> being the number of elements in x
 1226 
 1227 =item * the population RMS deviation from the mean:
 1228 
 1229   PRMS = sqrt( sum( (x-mean(x))^2 )/(N-1)
 1230 
 1231 The population deviation is the best-estimate of the deviation
 1232 of the population from which a sample is drawn.
 1233 
 1234 =item * the median
 1235 
 1236 The median is the 50th percentile data value.  Median is found by
 1237 L<medover|PDL::Ufunc/medover>, so WEIGHTING IS IGNORED FOR THE MEDIAN CALCULATION.
 1238 
 1239 =item * the minimum
 1240 
 1241 =item * the maximum
 1242 
 1243 =item * the average absolute deviation:
 1244 
 1245   AADEV = sum( abs(x-mean(x)) )/N
 1246 
 1247 =item * RMS deviation from the mean:
 1248 
 1249   RMS = sqrt(sum( (x-mean(x))^2 )/N)
 1250 
 1251 (also known as the root-mean-square deviation, or the square root of the
 1252 variance)
 1253 
 1254 =back
 1255 
 1256 This operator is a projection operator so the calculation
 1257 will take place over the final dimension. Thus if the input
 1258 is N-dimensional each returned value will be N-1 dimensional,
 1259 to calculate the statistics for the entire ndarray either
 1260 use C<clump(-1)> directly on the ndarray or call C<stats>.
 1261 
 1262 
 1263 
 1264 =for bad
 1265 
 1266 =for bad
 1267 
 1268 Bad values are simply ignored in the calculation, effectively reducing
 1269 the sample size.  If all data are bad then the output data are marked bad.
 1270 
 1271 
 1272 
 1273 =cut
 1274 #line 1275 "Primitive.pm"
 1275 
 1276 
 1277 
 1278 #line 949 "../../blib/lib/PDL/PP.pm"
 1279 
 1280 
 1281 sub PDL::statsover {
 1282    barf('Usage: ($mean,[$prms, $median, $min, $max, $adev, $rms]) = statsover($data,[$weights])') if @_>2;
 1283    my ($data, $weights) = @_;
 1284    $weights //= $data->ones();
 1285    my $median = $data->medover;
 1286    my $mean = PDL->nullcreate($data);
 1287    my $rms = PDL->nullcreate($data);
 1288    my $min = PDL->nullcreate($data);
 1289    my $max = PDL->nullcreate($data);
 1290    my $adev = PDL->nullcreate($data);
 1291    my $prms = PDL->nullcreate($data);
 1292    PDL::_statsover_int($data, $weights, $mean, $prms, $median, $min, $max, $adev, $rms);
 1293    wantarray ? ($mean, $prms, $median, $min, $max, $adev, $rms) : $mean;
 1294 }
 1295 #line 1296 "Primitive.pm"
 1296 
 1297 
 1298 
 1299 #line 950 "../../blib/lib/PDL/PP.pm"
 1300 
 1301 *statsover = \&PDL::statsover;
 1302 #line 1303 "Primitive.pm"
 1303 
 1304 
 1305 
 1306 #line 1184 "primitive.pd"
 1307 
 1308 =head2 stats
 1309 
 1310 =for ref
 1311 
 1312 Calculates useful statistics on an ndarray
 1313 
 1314 =for usage
 1315 
 1316  ($mean,$prms,$median,$min,$max,$adev,$rms) = stats($ndarray,[$weights]);
 1317 
 1318 This utility calculates all the most useful quantities in one call.
 1319 It works the same way as L</statsover>, except that the quantities are
 1320 calculated considering the entire input PDL as a single sample, rather
 1321 than as a collection of rows. See L</statsover> for definitions of the
 1322 returned quantities.
 1323 
 1324 =for bad
 1325 
 1326 Bad values are handled; if all input values are bad, then all of the output
 1327 values are flagged bad.
 1328 
 1329 =cut
 1330 
 1331 *stats    = \&PDL::stats;
 1332 sub PDL::stats {
 1333     barf('Usage: ($mean,[$rms]) = stats($data,[$weights])') if @_>2;
 1334     my ($data,$weights) = @_;
 1335 
 1336     # Ensure that $weights is properly broadcasted over; this could be
 1337     # done rather more efficiently...
 1338     if(defined $weights) {
 1339     $weights = pdl($weights) unless UNIVERSAL::isa($weights,'PDL');
 1340     if( ($weights->ndims != $data->ndims) or
 1341         (pdl($weights->dims) != pdl($data->dims))->or
 1342       ) {
 1343         $weights = $weights + zeroes($data)
 1344     }
 1345     $weights = $weights->flat;
 1346     }
 1347 
 1348     return PDL::statsover($data->flat,$weights);
 1349 }
 1350 #line 1351 "Primitive.pm"
 1351 
 1352 
 1353 
 1354 #line 948 "../../blib/lib/PDL/PP.pm"
 1355 
 1356 
 1357 
 1358 =head2 histogram
 1359 
 1360 =for sig
 1361 
 1362   Signature: (in(n); int+[o] hist(m); double step; double min; int msize => m)
 1363 
 1364 =for ref
 1365 
 1366 Calculates a histogram for given stepsize and minimum.
 1367 
 1368 =for usage
 1369 
 1370  $h = histogram($data, $step, $min, $numbins);
 1371  $hist = zeroes $numbins;  # Put histogram in existing ndarray.
 1372  histogram($data, $hist, $step, $min, $numbins);
 1373 
 1374 The histogram will contain C<$numbins> bins starting from C<$min>, each
 1375 C<$step> wide. The value in each bin is the number of
 1376 values in C<$data> that lie within the bin limits.
 1377 
 1378 
 1379 Data below the lower limit is put in the first bin, and data above the
 1380 upper limit is put in the last bin.
 1381 
 1382 The output is reset in a different broadcastloop so that you
 1383 can take a histogram of C<$a(10,12)> into C<$b(15)> and get the result
 1384 you want.
 1385 
 1386 For a higher-level interface, see L<hist|PDL::Basic/hist>.
 1387 
 1388 =for example
 1389 
 1390  pdl> p histogram(pdl(1,1,2),1,0,3)
 1391  [0 2 1]
 1392 
 1393 
 1394 
 1395 =for bad
 1396 
 1397 histogram processes bad values.
 1398 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1399 
 1400 
 1401 =cut
 1402 #line 1403 "Primitive.pm"
 1403 
 1404 
 1405 
 1406 #line 950 "../../blib/lib/PDL/PP.pm"
 1407 
 1408 *histogram = \&PDL::histogram;
 1409 #line 1410 "Primitive.pm"
 1410 
 1411 
 1412 
 1413 #line 948 "../../blib/lib/PDL/PP.pm"
 1414 
 1415 
 1416 
 1417 =head2 whistogram
 1418 
 1419 =for sig
 1420 
 1421   Signature: (in(n); float+ wt(n);float+[o] hist(m); double step; double min; int msize => m)
 1422 
 1423 =for ref
 1424 
 1425 Calculates a histogram from weighted data for given stepsize and minimum.
 1426 
 1427 =for usage
 1428 
 1429  $h = whistogram($data, $weights, $step, $min, $numbins);
 1430  $hist = zeroes $numbins;  # Put histogram in existing ndarray.
 1431  whistogram($data, $weights, $hist, $step, $min, $numbins);
 1432 
 1433 The histogram will contain C<$numbins> bins starting from C<$min>, each
 1434 C<$step> wide. The value in each bin is the sum of the values in C<$weights>
 1435 that correspond to values in C<$data> that lie within the bin limits.
 1436 
 1437 Data below the lower limit is put in the first bin, and data above the
 1438 upper limit is put in the last bin.
 1439 
 1440 The output is reset in a different broadcastloop so that you
 1441 can take a histogram of C<$a(10,12)> into C<$b(15)> and get the result
 1442 you want.
 1443 
 1444 =for example
 1445 
 1446  pdl> p whistogram(pdl(1,1,2), pdl(0.1,0.1,0.5), 1, 0, 4)
 1447  [0 0.2 0.5 0]
 1448 
 1449 
 1450 
 1451 =for bad
 1452 
 1453 whistogram processes bad values.
 1454 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1455 
 1456 
 1457 =cut
 1458 #line 1459 "Primitive.pm"
 1459 
 1460 
 1461 
 1462 #line 950 "../../blib/lib/PDL/PP.pm"
 1463 
 1464 *whistogram = \&PDL::whistogram;
 1465 #line 1466 "Primitive.pm"
 1466 
 1467 
 1468 
 1469 #line 948 "../../blib/lib/PDL/PP.pm"
 1470 
 1471 
 1472 
 1473 =head2 histogram2d
 1474 
 1475 =for sig
 1476 
 1477   Signature: (ina(n); inb(n); int+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
 1478                  double stepb; double minb; int mbsize => mb;)
 1479 
 1480 =for ref
 1481 
 1482 Calculates a 2d histogram.
 1483 
 1484 =for usage
 1485 
 1486  $h = histogram2d($datax, $datay, $stepx, $minx,
 1487        $nbinx, $stepy, $miny, $nbiny);
 1488  $hist = zeroes $nbinx, $nbiny;  # Put histogram in existing ndarray.
 1489  histogram2d($datax, $datay, $hist, $stepx, $minx,
 1490        $nbinx, $stepy, $miny, $nbiny);
 1491 
 1492 The histogram will contain C<$nbinx> x C<$nbiny> bins, with the lower
 1493 limits of the first one at C<($minx, $miny)>, and with bin size
 1494 C<($stepx, $stepy)>.
 1495 The value in each bin is the number of
 1496 values in C<$datax> and C<$datay> that lie within the bin limits.
 1497 
 1498 Data below the lower limit is put in the first bin, and data above the
 1499 upper limit is put in the last bin.
 1500 
 1501 =for example
 1502 
 1503  pdl> p histogram2d(pdl(1,1,1,2,2),pdl(2,1,1,1,1),1,0,3,1,0,3)
 1504  [
 1505   [0 0 0]
 1506   [0 2 2]
 1507   [0 1 0]
 1508  ]
 1509 
 1510 
 1511 
 1512 =for bad
 1513 
 1514 histogram2d processes bad values.
 1515 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1516 
 1517 
 1518 =cut
 1519 #line 1520 "Primitive.pm"
 1520 
 1521 
 1522 
 1523 #line 950 "../../blib/lib/PDL/PP.pm"
 1524 
 1525 *histogram2d = \&PDL::histogram2d;
 1526 #line 1527 "Primitive.pm"
 1527 
 1528 
 1529 
 1530 #line 948 "../../blib/lib/PDL/PP.pm"
 1531 
 1532 
 1533 
 1534 =head2 whistogram2d
 1535 
 1536 =for sig
 1537 
 1538   Signature: (ina(n); inb(n); float+ wt(n);float+[o] hist(ma,mb); double stepa; double mina; int masize => ma;
 1539                  double stepb; double minb; int mbsize => mb;)
 1540 
 1541 =for ref
 1542 
 1543 Calculates a 2d histogram from weighted data.
 1544 
 1545 =for usage
 1546 
 1547  $h = whistogram2d($datax, $datay, $weights,
 1548        $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
 1549  $hist = zeroes $nbinx, $nbiny;  # Put histogram in existing ndarray.
 1550  whistogram2d($datax, $datay, $weights, $hist,
 1551        $stepx, $minx, $nbinx, $stepy, $miny, $nbiny);
 1552 
 1553 The histogram will contain C<$nbinx> x C<$nbiny> bins, with the lower
 1554 limits of the first one at C<($minx, $miny)>, and with bin size
 1555 C<($stepx, $stepy)>.
 1556 The value in each bin is the sum of the values in
 1557 C<$weights> that correspond to values in C<$datax> and C<$datay> that lie within the bin limits.
 1558 
 1559 Data below the lower limit is put in the first bin, and data above the
 1560 upper limit is put in the last bin.
 1561 
 1562 =for example
 1563 
 1564  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)
 1565  [
 1566   [  0   0   0]
 1567   [  0 0.5 0.9]
 1568   [  0 0.1   0]
 1569  ]
 1570 
 1571 
 1572 
 1573 =for bad
 1574 
 1575 whistogram2d processes bad values.
 1576 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1577 
 1578 
 1579 =cut
 1580 #line 1581 "Primitive.pm"
 1581 
 1582 
 1583 
 1584 #line 950 "../../blib/lib/PDL/PP.pm"
 1585 
 1586 *whistogram2d = \&PDL::whistogram2d;
 1587 #line 1588 "Primitive.pm"
 1588 
 1589 
 1590 
 1591 #line 948 "../../blib/lib/PDL/PP.pm"
 1592 
 1593 
 1594 
 1595 =head2 fibonacci
 1596 
 1597 =for sig
 1598 
 1599   Signature: (i(n); indx [o]x(n))
 1600 
 1601 =for ref
 1602 
 1603 Constructor - a vector with Fibonacci's sequence
 1604 
 1605 =for bad
 1606 
 1607 fibonacci does not process bad values.
 1608 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1609 
 1610 
 1611 =cut
 1612 #line 1613 "Primitive.pm"
 1613 
 1614 
 1615 
 1616 #line 949 "../../blib/lib/PDL/PP.pm"
 1617 
 1618 sub fibonacci { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->fibonacci : PDL->fibonacci(@_) }
 1619 sub PDL::fibonacci{
 1620    my $x = &PDL::Core::_construct;
 1621    my $is_inplace = $x->is_inplace;
 1622    my ($in, $out) = $x->clump(-1);
 1623    $out = $is_inplace ? $in->inplace : PDL->null;
 1624    PDL::_fibonacci_int($in, $out);
 1625    $out;
 1626 }
 1627 #line 1628 "Primitive.pm"
 1628 
 1629 
 1630 
 1631 #line 950 "../../blib/lib/PDL/PP.pm"
 1632 #line 1633 "Primitive.pm"
 1633 
 1634 
 1635 
 1636 #line 948 "../../blib/lib/PDL/PP.pm"
 1637 
 1638 
 1639 
 1640 =head2 append
 1641 
 1642 =for sig
 1643 
 1644   Signature: (a(n); b(m); [o] c(mn))
 1645 
 1646 
 1647 =for ref
 1648 
 1649 append two ndarrays by concatenating along their first dimensions
 1650 
 1651 =for example
 1652 
 1653  $x = ones(2,4,7);
 1654  $y = sequence 5;
 1655  $c = $x->append($y);  # size of $c is now (7,4,7) (a jumbo-ndarray ;)
 1656 
 1657 C<append> appends two ndarrays along their first dimensions. The rest of the
 1658 dimensions must be compatible in the broadcasting sense. The resulting
 1659 size of the first dimension is the sum of the sizes of the first dimensions
 1660 of the two argument ndarrays - i.e. C<n + m>.
 1661 
 1662 Similar functions include L</glue> (below), which can append more
 1663 than two ndarrays along an arbitrary dimension, and
 1664 L<cat|PDL::Core/cat>, which can append more than two ndarrays that all
 1665 have the same sized dimensions.
 1666 
 1667 
 1668 =for bad
 1669 
 1670 append does not process bad values.
 1671 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1672 
 1673 
 1674 =cut
 1675 #line 1676 "Primitive.pm"
 1676 
 1677 
 1678 
 1679 #line 949 "../../blib/lib/PDL/PP.pm"
 1680 
 1681 
 1682 #line 1491 "primitive.pd"
 1683 
 1684 sub PDL::append {
 1685   my ($i1, $i2, $o) = map PDL->topdl($_), @_;
 1686   if (grep $_->isempty, $i1, $i2) {
 1687     if (!defined $o) {
 1688       return $i2->isnull ? PDL->zeroes(0) : $i2->copy if $i1->isempty;
 1689       return $i1->isnull ? PDL->zeroes(0) : $i1->copy;
 1690     } else {
 1691       $o .= $i2->isnull ? PDL->zeroes(0) : $i2, return $o if $i1->isempty;
 1692       $o .= $i1->isnull ? PDL->zeroes(0) : $i1, return $o;
 1693     }
 1694   }
 1695   $o //= PDL->null;
 1696   PDL::_append_int($i1, $i2->convert($i1->type), $o);
 1697   $o;
 1698 }
 1699         
 1700 #line 969 "../../blib/lib/PDL/PP.pm"
 1701 #line 1702 "Primitive.pm"
 1702 
 1703 
 1704 
 1705 #line 950 "../../blib/lib/PDL/PP.pm"
 1706 
 1707 *append = \&PDL::append;
 1708 #line 1709 "Primitive.pm"
 1709 
 1710 
 1711 
 1712 #line 1538 "primitive.pd"
 1713 
 1714 =head2 glue
 1715 
 1716 =for usage
 1717 
 1718   $c = $x->glue(<dim>,$y,...)
 1719 
 1720 =for ref
 1721 
 1722 Glue two or more PDLs together along an arbitrary dimension
 1723 (N-D L</append>).
 1724 
 1725 Sticks $x, $y, and all following arguments together along the
 1726 specified dimension.  All other dimensions must be compatible in the
 1727 broadcasting sense.
 1728 
 1729 Glue is permissive, in the sense that every PDL is treated as having an
 1730 infinite number of trivial dimensions of order 1 -- so C<< $x->glue(3,$y) >>
 1731 works, even if $x and $y are only one dimensional.
 1732 
 1733 If one of the PDLs has no elements, it is ignored.  Likewise, if one
 1734 of them is actually the undefined value, it is treated as if it had no
 1735 elements.
 1736 
 1737 If the first parameter is a defined perl scalar rather than a pdl,
 1738 then it is taken as a dimension along which to glue everything else,
 1739 so you can say C<$cube = PDL::glue(3,@image_list);> if you like.
 1740 
 1741 C<glue> is implemented in pdl, using a combination of L<xchg|PDL::Slices/xchg> and
 1742 L</append>.  It should probably be updated (one day) to a pure PP
 1743 function.
 1744 
 1745 Similar functions include L</append> (above), which appends
 1746 only two ndarrays along their first dimension, and
 1747 L<cat|PDL::Core/cat>, which can append more than two ndarrays that all
 1748 have the same sized dimensions.
 1749 
 1750 =cut
 1751 
 1752 sub PDL::glue{
 1753     my($x) = shift;
 1754     my($dim) = shift;
 1755 
 1756     ($dim, $x) = ($x, $dim) if defined $x && !ref $x;
 1757     confess 'dimension must be Perl scalar' if ref $dim;
 1758 
 1759     if(!defined $x || $x->nelem==0) {
 1760     return $x unless(@_);
 1761     return shift() if(@_<=1);
 1762     $x=shift;
 1763     return PDL::glue($x,$dim,@_);
 1764     }
 1765 
 1766     if($dim - $x->dim(0) > 100) {
 1767     print STDERR "warning:: PDL::glue allocating >100 dimensions!\n";
 1768     }
 1769     while($dim >= $x->ndims) {
 1770     $x = $x->dummy(-1,1);
 1771     }
 1772     $x = $x->xchg(0,$dim);
 1773 
 1774     while(scalar(@_)){
 1775     my $y = shift;
 1776     next unless(defined $y && $y->nelem);
 1777 
 1778     while($dim >= $y->ndims) {
 1779         $y = $y->dummy(-1,1);
 1780         }
 1781     $y = $y->xchg(0,$dim);
 1782     $x = $x->append($y);
 1783     }
 1784     $x->xchg(0,$dim);
 1785 }
 1786 #line 1787 "Primitive.pm"
 1787 
 1788 
 1789 
 1790 #line 950 "../../blib/lib/PDL/PP.pm"
 1791 
 1792 *axisvalues = \&PDL::axisvalues;
 1793 #line 1794 "Primitive.pm"
 1794 
 1795 
 1796 
 1797 #line 948 "../../blib/lib/PDL/PP.pm"
 1798 
 1799 
 1800 
 1801 =head2 cmpvec
 1802 
 1803 =for sig
 1804 
 1805   Signature: (a(n); b(n); sbyte [o]c())
 1806 
 1807 
 1808 =for ref
 1809 
 1810 Compare two vectors lexicographically.
 1811 
 1812 Returns -1 if a is less, 1 if greater, 0 if equal.
 1813 
 1814 
 1815 =for bad
 1816 
 1817 The output is bad if any input values up to the point of inequality are
 1818 bad - any after are ignored.
 1819 
 1820 
 1821 =cut
 1822 #line 1823 "Primitive.pm"
 1823 
 1824 
 1825 
 1826 #line 950 "../../blib/lib/PDL/PP.pm"
 1827 
 1828 *cmpvec = \&PDL::cmpvec;
 1829 #line 1830 "Primitive.pm"
 1830 
 1831 
 1832 
 1833 #line 948 "../../blib/lib/PDL/PP.pm"
 1834 
 1835 
 1836 
 1837 =head2 eqvec
 1838 
 1839 =for sig
 1840 
 1841   Signature: (a(n); b(n); sbyte [o]c())
 1842 
 1843 
 1844 =for ref
 1845 
 1846 Compare two vectors, returning 1 if equal, 0 if not equal.
 1847 
 1848 
 1849 =for bad
 1850 
 1851 The output is bad if any input values are bad.
 1852 
 1853 =cut
 1854 #line 1855 "Primitive.pm"
 1855 
 1856 
 1857 
 1858 #line 950 "../../blib/lib/PDL/PP.pm"
 1859 
 1860 *eqvec = \&PDL::eqvec;
 1861 #line 1862 "Primitive.pm"
 1862 
 1863 
 1864 
 1865 #line 948 "../../blib/lib/PDL/PP.pm"
 1866 
 1867 
 1868 
 1869 =head2 enumvec
 1870 
 1871 =for sig
 1872 
 1873   Signature: (v(M,N); indx [o]k(N))
 1874 
 1875 =for ref
 1876 
 1877 Enumerate a list of vectors with locally unique keys.
 1878 
 1879 Given a sorted list of vectors $v, generate a vector $k containing locally unique keys for the elements of $v
 1880 (where an "element" is a vector of length $M ocurring in $v).
 1881 
 1882 Note that the keys returned in $k are only unique over a run of a single vector in $v,
 1883 so that each unique vector in $v has at least one 0 (zero) index in $k associated with it.
 1884 If you need global keys, see enumvecg().
 1885 
 1886 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1887 
 1888 
 1889 =for bad
 1890 
 1891 enumvec does not process bad values.
 1892 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1893 
 1894 
 1895 =cut
 1896 #line 1897 "Primitive.pm"
 1897 
 1898 
 1899 
 1900 #line 950 "../../blib/lib/PDL/PP.pm"
 1901 
 1902 *enumvec = \&PDL::enumvec;
 1903 #line 1904 "Primitive.pm"
 1904 
 1905 
 1906 
 1907 #line 948 "../../blib/lib/PDL/PP.pm"
 1908 
 1909 
 1910 
 1911 =head2 enumvecg
 1912 
 1913 =for sig
 1914 
 1915   Signature: (v(M,N); indx [o]k(N))
 1916 
 1917 =for ref
 1918 
 1919 Enumerate a list of vectors with globally unique keys.
 1920 
 1921 Given a sorted list of vectors $v, generate a vector $k containing globally unique keys for the elements of $v
 1922 (where an "element" is a vector of length $M ocurring in $v).
 1923 Basically does the same thing as:
 1924 
 1925  $k = $v->vsearchvec($v->uniqvec);
 1926 
 1927 ... but somewhat more efficiently.
 1928 
 1929 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1930 
 1931 
 1932 =for bad
 1933 
 1934 enumvecg does not process bad values.
 1935 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1936 
 1937 
 1938 =cut
 1939 #line 1940 "Primitive.pm"
 1940 
 1941 
 1942 
 1943 #line 950 "../../blib/lib/PDL/PP.pm"
 1944 
 1945 *enumvecg = \&PDL::enumvecg;
 1946 #line 1947 "Primitive.pm"
 1947 
 1948 
 1949 
 1950 #line 948 "../../blib/lib/PDL/PP.pm"
 1951 
 1952 
 1953 
 1954 =head2 vsearchvec
 1955 
 1956 =for sig
 1957 
 1958   Signature: (find(M); which(M,N); indx [o]found())
 1959 
 1960 =for ref
 1961 
 1962 Routine for searching N-dimensional values - akin to vsearch() for vectors.
 1963 
 1964 =for usage
 1965 
 1966  $found   = vsearchvec($find, $which);
 1967  $nearest = $which->dice_axis(1,$found);
 1968 
 1969 Returns for each row-vector in C<$find> the index along dimension N
 1970 of the least row vector of C<$which>
 1971 greater or equal to it.
 1972 C<$which> should be sorted in increasing order.
 1973 If the value of C<$find> is larger
 1974 than any member of C<$which>, the index to the last element of C<$which> is
 1975 returned.
 1976 
 1977 See also: L</vsearch>.
 1978 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 1979 
 1980 
 1981 =for bad
 1982 
 1983 vsearchvec does not process bad values.
 1984 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 1985 
 1986 
 1987 =cut
 1988 #line 1989 "Primitive.pm"
 1989 
 1990 
 1991 
 1992 #line 950 "../../blib/lib/PDL/PP.pm"
 1993 
 1994 *vsearchvec = \&PDL::vsearchvec;
 1995 #line 1996 "Primitive.pm"
 1996 
 1997 
 1998 
 1999 #line 948 "../../blib/lib/PDL/PP.pm"
 2000 
 2001 
 2002 
 2003 =head2 unionvec
 2004 
 2005 =for sig
 2006 
 2007   Signature: (a(M,NA); b(M,NB); [o]c(M,NC); indx [o]nc())
 2008 
 2009 =for ref
 2010 
 2011 Union of two vector-valued PDLs.
 2012 
 2013 Input PDLs $a() and $b() B<MUST> be sorted in lexicographic order.
 2014 On return, $nc() holds the actual number of vector-values in the union.
 2015 
 2016 In scalar context, slices $c() to the actual number of elements in the union
 2017 and returns the sliced PDL.
 2018 
 2019 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 2020 
 2021 
 2022 =for bad
 2023 
 2024 unionvec does not process bad values.
 2025 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2026 
 2027 
 2028 =cut
 2029 #line 2030 "Primitive.pm"
 2030 
 2031 
 2032 
 2033 #line 949 "../../blib/lib/PDL/PP.pm"
 2034 
 2035 
 2036  sub PDL::unionvec {
 2037    my ($a,$b,$c,$nc) = @_;
 2038    $c = PDL->null if (!defined($nc));
 2039    $nc = PDL->null if (!defined($nc));
 2040    PDL::_unionvec_int($a,$b,$c,$nc);
 2041    return ($c,$nc) if (wantarray);
 2042    return $c->slice(",0:".($nc->max-1));
 2043  }
 2044 #line 2045 "Primitive.pm"
 2045 
 2046 
 2047 
 2048 #line 950 "../../blib/lib/PDL/PP.pm"
 2049 
 2050 *unionvec = \&PDL::unionvec;
 2051 #line 2052 "Primitive.pm"
 2052 
 2053 
 2054 
 2055 #line 948 "../../blib/lib/PDL/PP.pm"
 2056 
 2057 
 2058 
 2059 =head2 intersectvec
 2060 
 2061 =for sig
 2062 
 2063   Signature: (a(M,NA); b(M,NB); [o]c(M,NC); indx [o]nc())
 2064 
 2065 =for ref
 2066 
 2067 Intersection of two vector-valued PDLs.
 2068 Input PDLs $a() and $b() B<MUST> be sorted in lexicographic order.
 2069 On return, $nc() holds the actual number of vector-values in the intersection.
 2070 
 2071 In scalar context, slices $c() to the actual number of elements in the intersection
 2072 and returns the sliced PDL.
 2073 
 2074 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 2075 
 2076 
 2077 =for bad
 2078 
 2079 intersectvec does not process bad values.
 2080 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2081 
 2082 
 2083 =cut
 2084 #line 2085 "Primitive.pm"
 2085 
 2086 
 2087 
 2088 #line 949 "../../blib/lib/PDL/PP.pm"
 2089 
 2090 
 2091  sub PDL::intersectvec {
 2092    my ($a,$b,$c,$nc) = @_;
 2093    $c = PDL->null if (!defined($c));
 2094    $nc = PDL->null if (!defined($nc));
 2095    PDL::_intersectvec_int($a,$b,$c,$nc);
 2096    return ($c,$nc) if (wantarray);
 2097    my $nc_max = $nc->max;
 2098    return ($nc_max > 0
 2099        ? $c->slice(",0:".($nc_max-1))
 2100        : $c->reshape($c->dim(0), 0, ($c->dims)[2..($c->ndims-1)]));
 2101  }
 2102 #line 2103 "Primitive.pm"
 2103 
 2104 
 2105 
 2106 #line 950 "../../blib/lib/PDL/PP.pm"
 2107 
 2108 *intersectvec = \&PDL::intersectvec;
 2109 #line 2110 "Primitive.pm"
 2110 
 2111 
 2112 
 2113 #line 948 "../../blib/lib/PDL/PP.pm"
 2114 
 2115 
 2116 
 2117 =head2 setdiffvec
 2118 
 2119 =for sig
 2120 
 2121   Signature: (a(M,NA); b(M,NB); [o]c(M,NC); indx [o]nc())
 2122 
 2123 =for ref
 2124 
 2125 Set-difference ($a() \ $b()) of two vector-valued PDLs.
 2126 
 2127 Input PDLs $a() and $b() B<MUST> be sorted in lexicographic order.
 2128 On return, $nc() holds the actual number of vector-values in the computed vector set.
 2129 
 2130 In scalar context, slices $c() to the actual number of elements in the output vector set
 2131 and returns the sliced PDL.
 2132 
 2133 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 2134 
 2135 
 2136 =for bad
 2137 
 2138 setdiffvec does not process bad values.
 2139 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2140 
 2141 
 2142 =cut
 2143 #line 2144 "Primitive.pm"
 2144 
 2145 
 2146 
 2147 #line 949 "../../blib/lib/PDL/PP.pm"
 2148 
 2149 
 2150  sub PDL::setdiffvec {
 2151   my ($a,$b,$c,$nc) = @_;
 2152   $c = PDL->null if (!defined($c));
 2153   $nc = PDL->null if (!defined($nc));
 2154   PDL::_setdiffvec_int($a,$b,$c,$nc);
 2155   return ($c,$nc) if (wantarray);
 2156   my $nc_max = $nc->max;
 2157   return ($nc_max > 0
 2158       ? $c->slice(",0:".($nc_max-1))
 2159       : $c->reshape($c->dim(0), 0, ($c->dims)[2..($c->ndims-1)]));
 2160  }
 2161 #line 2162 "Primitive.pm"
 2162 
 2163 
 2164 
 2165 #line 950 "../../blib/lib/PDL/PP.pm"
 2166 
 2167 *setdiffvec = \&PDL::setdiffvec;
 2168 #line 2169 "Primitive.pm"
 2169 
 2170 
 2171 
 2172 #line 948 "../../blib/lib/PDL/PP.pm"
 2173 
 2174 
 2175 
 2176 =head2 union_sorted
 2177 
 2178 =for sig
 2179 
 2180   Signature: (a(NA); b(NB); [o]c(NC); indx [o]nc())
 2181 
 2182 =for ref
 2183 
 2184 Union of two flat sorted unique-valued PDLs.
 2185 Input PDLs $a() and $b() B<MUST> be sorted in lexicographic order and contain no duplicates.
 2186 On return, $nc() holds the actual number of values in the union.
 2187 
 2188 In scalar context, reshapes $c() to the actual number of elements in the union and returns it.
 2189 
 2190 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 2191 
 2192 
 2193 =for bad
 2194 
 2195 union_sorted does not process bad values.
 2196 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2197 
 2198 
 2199 =cut
 2200 #line 2201 "Primitive.pm"
 2201 
 2202 
 2203 
 2204 #line 949 "../../blib/lib/PDL/PP.pm"
 2205 
 2206 
 2207  sub PDL::union_sorted {
 2208    my ($a,$b,$c,$nc) = @_;
 2209    $c = PDL->null if (!defined($c));
 2210    $nc = PDL->null if (!defined($nc));
 2211    PDL::_union_sorted_int($a,$b,$c,$nc);
 2212    return ($c,$nc) if (wantarray);
 2213    return $c->slice("0:".($nc->max-1));
 2214  }
 2215 #line 2216 "Primitive.pm"
 2216 
 2217 
 2218 
 2219 #line 950 "../../blib/lib/PDL/PP.pm"
 2220 
 2221 *union_sorted = \&PDL::union_sorted;
 2222 #line 2223 "Primitive.pm"
 2223 
 2224 
 2225 
 2226 #line 948 "../../blib/lib/PDL/PP.pm"
 2227 
 2228 
 2229 
 2230 =head2 intersect_sorted
 2231 
 2232 =for sig
 2233 
 2234   Signature: (a(NA); b(NB); [o]c(NC); indx [o]nc())
 2235 
 2236 =for ref
 2237 
 2238 Intersection of two flat sorted unique-valued PDLs.
 2239 Input PDLs $a() and $b() B<MUST> be sorted in lexicographic order and contain no duplicates.
 2240 On return, $nc() holds the actual number of values in the intersection.
 2241 
 2242 In scalar context, reshapes $c() to the actual number of elements in the intersection and returns it.
 2243 
 2244 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 2245 
 2246 
 2247 =for bad
 2248 
 2249 intersect_sorted does not process bad values.
 2250 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2251 
 2252 
 2253 =cut
 2254 #line 2255 "Primitive.pm"
 2255 
 2256 
 2257 
 2258 #line 949 "../../blib/lib/PDL/PP.pm"
 2259 
 2260 
 2261  sub PDL::intersect_sorted {
 2262    my ($a,$b,$c,$nc) = @_;
 2263    $c = PDL->null if (!defined($c));
 2264    $nc = PDL->null if (!defined($nc));
 2265    PDL::_intersect_sorted_int($a,$b,$c,$nc);
 2266    return ($c,$nc) if (wantarray);
 2267    my $nc_max = $nc->max;
 2268    return ($nc_max > 0
 2269        ? $c->slice("0:".($nc_max-1))
 2270        : $c->reshape(0, ($c->dims)[1..($c->ndims-1)]));
 2271  }
 2272 #line 2273 "Primitive.pm"
 2273 
 2274 
 2275 
 2276 #line 950 "../../blib/lib/PDL/PP.pm"
 2277 
 2278 *intersect_sorted = \&PDL::intersect_sorted;
 2279 #line 2280 "Primitive.pm"
 2280 
 2281 
 2282 
 2283 #line 948 "../../blib/lib/PDL/PP.pm"
 2284 
 2285 
 2286 
 2287 =head2 setdiff_sorted
 2288 
 2289 =for sig
 2290 
 2291   Signature: (a(NA); b(NB); [o]c(NC); indx [o]nc())
 2292 
 2293 =for ref
 2294 
 2295 Set-difference ($a() \ $b()) of two flat sorted unique-valued PDLs.
 2296 
 2297 Input PDLs $a() and $b() B<MUST> be sorted in lexicographic order and contain no duplicate values.
 2298 On return, $nc() holds the actual number of values in the computed vector set.
 2299 
 2300 In scalar context, reshapes $c() to the actual number of elements in the difference set and returns it.
 2301 
 2302 Contributed by Bryan Jurish E<lt>moocow@cpan.orgE<gt>.
 2303 
 2304 
 2305 =for bad
 2306 
 2307 setdiff_sorted does not process bad values.
 2308 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2309 
 2310 
 2311 =cut
 2312 #line 2313 "Primitive.pm"
 2313 
 2314 
 2315 
 2316 #line 949 "../../blib/lib/PDL/PP.pm"
 2317 
 2318 
 2319  sub PDL::setdiff_sorted {
 2320    my ($a,$b,$c,$nc) = @_;
 2321    $c = PDL->null if (!defined($c));
 2322    $nc = PDL->null if (!defined($nc));
 2323    PDL::_setdiff_sorted_int($a,$b,$c,$nc);
 2324    return ($c,$nc) if (wantarray);
 2325    my $nc_max = $nc->max;
 2326    return ($nc_max > 0
 2327        ? $c->slice("0:".($nc_max-1))
 2328        : $c->reshape(0, ($c->dims)[1..($c->ndims-1)]));
 2329  }
 2330 #line 2331 "Primitive.pm"
 2331 
 2332 
 2333 
 2334 #line 950 "../../blib/lib/PDL/PP.pm"
 2335 
 2336 *setdiff_sorted = \&PDL::setdiff_sorted;
 2337 #line 2338 "Primitive.pm"
 2338 
 2339 
 2340 
 2341 #line 948 "../../blib/lib/PDL/PP.pm"
 2342 
 2343 
 2344 
 2345 =head2 srand
 2346 
 2347 =for sig
 2348 
 2349   Signature: (a())
 2350 
 2351 =for ref
 2352 
 2353 Seed random-number generator with a 64-bit int. Will generate seed data
 2354 for a number of threads equal to the return-value of
 2355 L<PDL::Core/online_cpus>.
 2356 
 2357 =for usage
 2358 
 2359  srand(); # uses current time
 2360  srand(5); # fixed number e.g. for testing
 2361 
 2362 
 2363 
 2364 =for bad
 2365 
 2366 srand does not process bad values.
 2367 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2368 
 2369 
 2370 =cut
 2371 #line 2372 "Primitive.pm"
 2372 
 2373 
 2374 
 2375 #line 949 "../../blib/lib/PDL/PP.pm"
 2376 
 2377 *srand = \&PDL::srand;
 2378 sub PDL::srand { PDL::_srand_int($_[0] // PDL::Core::seed()) }
 2379 #line 2380 "Primitive.pm"
 2380 
 2381 
 2382 
 2383 #line 950 "../../blib/lib/PDL/PP.pm"
 2384 
 2385 *srand = \&PDL::srand;
 2386 #line 2387 "Primitive.pm"
 2387 
 2388 
 2389 
 2390 #line 948 "../../blib/lib/PDL/PP.pm"
 2391 
 2392 
 2393 
 2394 =head2 random
 2395 
 2396 =for sig
 2397 
 2398   Signature: (a())
 2399 
 2400 =for ref
 2401 
 2402 Constructor which returns ndarray of random numbers
 2403 
 2404 =for usage
 2405 
 2406  $x = random([type], $nx, $ny, $nz,...);
 2407  $x = random $y;
 2408 
 2409 etc (see L<zeroes|PDL::Core/zeroes>).
 2410 
 2411 This is the uniform distribution between 0 and 1 (assumedly
 2412 excluding 1 itself). The arguments are the same as C<zeroes>
 2413 (q.v.) - i.e. one can specify dimensions, types or give
 2414 a template.
 2415 
 2416 You can use the PDL function L</srand> to seed the random generator.
 2417 If it has not been called yet, it will be with the current time.
 2418 
 2419 
 2420 =for bad
 2421 
 2422 random does not process bad values.
 2423 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2424 
 2425 
 2426 =cut
 2427 #line 2428 "Primitive.pm"
 2428 
 2429 
 2430 
 2431 #line 949 "../../blib/lib/PDL/PP.pm"
 2432 
 2433 sub random { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->random : PDL->random(@_) }
 2434 sub PDL::random {
 2435    my $class = shift;
 2436    my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
 2437    PDL::_random_int($x);
 2438    return $x;
 2439 }
 2440 #line 2441 "Primitive.pm"
 2441 
 2442 
 2443 
 2444 #line 950 "../../blib/lib/PDL/PP.pm"
 2445 #line 2446 "Primitive.pm"
 2446 
 2447 
 2448 
 2449 #line 948 "../../blib/lib/PDL/PP.pm"
 2450 
 2451 
 2452 
 2453 =head2 randsym
 2454 
 2455 =for sig
 2456 
 2457   Signature: (a())
 2458 
 2459 =for ref
 2460 
 2461 Constructor which returns ndarray of random numbers
 2462 
 2463 =for usage
 2464 
 2465  $x = randsym([type], $nx, $ny, $nz,...);
 2466  $x = randsym $y;
 2467 
 2468 etc (see L<zeroes|PDL::Core/zeroes>).
 2469 
 2470 This is the uniform distribution between 0 and 1 (excluding both 0 and
 2471 1, cf L</random>). The arguments are the same as C<zeroes> (q.v.) -
 2472 i.e. one can specify dimensions, types or give a template.
 2473 
 2474 You can use the PDL function L</srand> to seed the random generator.
 2475 If it has not been called yet, it will be with the current time.
 2476 
 2477 
 2478 =for bad
 2479 
 2480 randsym does not process bad values.
 2481 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 2482 
 2483 
 2484 =cut
 2485 #line 2486 "Primitive.pm"
 2486 
 2487 
 2488 
 2489 #line 949 "../../blib/lib/PDL/PP.pm"
 2490 
 2491 sub randsym { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->randsym : PDL->randsym(@_) }
 2492 sub PDL::randsym {
 2493    my $class = shift;
 2494    my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
 2495    PDL::_randsym_int($x);
 2496    return $x;
 2497 }
 2498 #line 2499 "Primitive.pm"
 2499 
 2500 
 2501 
 2502 #line 950 "../../blib/lib/PDL/PP.pm"
 2503 #line 2504 "Primitive.pm"
 2504 
 2505 
 2506 
 2507 #line 2313 "primitive.pd"
 2508 
 2509 =head2 grandom
 2510 
 2511 =for ref
 2512 
 2513 Constructor which returns ndarray of Gaussian random numbers
 2514 
 2515 =for usage
 2516 
 2517  $x = grandom([type], $nx, $ny, $nz,...);
 2518  $x = grandom $y;
 2519 
 2520 etc (see L<zeroes|PDL::Core/zeroes>).
 2521 
 2522 This is generated using the math library routine C<ndtri>.
 2523 
 2524 Mean = 0, Stddev = 1
 2525 
 2526 You can use the PDL function L</srand> to seed the random generator.
 2527 If it has not been called yet, it will be with the current time.
 2528 
 2529 =cut
 2530 
 2531 sub grandom { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? $_[0]->grandom : PDL->grandom(@_) }
 2532 sub PDL::grandom {
 2533    my $class = shift;
 2534    my $x = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
 2535    use PDL::Math 'ndtri';
 2536    $x .= ndtri(randsym($x));
 2537    return $x;
 2538 }
 2539 #line 2540 "Primitive.pm"
 2540 
 2541 
 2542 
 2543 #line 2353 "primitive.pd"
 2544 
 2545 =head2 vsearch
 2546 
 2547 =for sig
 2548 
 2549   Signature: ( vals(); xs(n); [o] indx(); [\%options] )
 2550 
 2551 =for ref
 2552 
 2553 Efficiently search for values in a sorted ndarray, returning indices.
 2554 
 2555 =for usage
 2556 
 2557   $idx = vsearch( $vals, $x, [\%options] );
 2558   vsearch( $vals, $x, $idx, [\%options ] );
 2559 
 2560 B<vsearch> performs a binary search in the ordered ndarray C<$x>,
 2561 for the values from C<$vals> ndarray, returning indices into C<$x>.
 2562 What is a "match", and the meaning of the returned indices, are determined
 2563 by the options.
 2564 
 2565 The C<mode> option indicates which method of searching to use, and may
 2566 be one of:
 2567 
 2568 =over
 2569 
 2570 =item C<sample>
 2571 
 2572 invoke L<B<vsearch_sample>|/vsearch_sample>, returning indices appropriate for sampling
 2573 within a distribution.
 2574 
 2575 =item C<insert_leftmost>
 2576 
 2577 invoke L<B<vsearch_insert_leftmost>|/vsearch_insert_leftmost>, returning the left-most possible
 2578 insertion point which still leaves the ndarray sorted.
 2579 
 2580 =item C<insert_rightmost>
 2581 
 2582 invoke L<B<vsearch_insert_rightmost>|/vsearch_insert_rightmost>, returning the right-most possible
 2583 insertion point which still leaves the ndarray sorted.
 2584 
 2585 =item C<match>
 2586 
 2587 invoke L<B<vsearch_match>|/vsearch_match>, returning the index of a matching element,
 2588 else -(insertion point + 1)
 2589 
 2590 =item C<bin_inclusive>
 2591 
 2592 invoke L<B<vsearch_bin_inclusive>|/vsearch_bin_inclusive>, returning an index appropriate for binning
 2593 on a grid where the left bin edges are I<inclusive> of the bin. See
 2594 below for further explanation of the bin.
 2595 
 2596 =item C<bin_exclusive>
 2597 
 2598 invoke L<B<vsearch_bin_exclusive>|/vsearch_bin_exclusive>, returning an index appropriate for binning
 2599 on a grid where the left bin edges are I<exclusive> of the bin. See
 2600 below for further explanation of the bin.
 2601 
 2602 =back
 2603 
 2604 The default value of C<mode> is C<sample>.
 2605 
 2606 =for example
 2607 
 2608   use PDL;
 2609   
 2610   my @modes = qw( sample insert_leftmost insert_rightmost match
 2611                   bin_inclusive bin_exclusive );
 2612   
 2613   # Generate a sequence of 3 zeros, 3 ones, ..., 3 fours.
 2614   my $x = zeroes(3,5)->yvals->flat;
 2615   
 2616   for my $mode ( @modes ) {
 2617     # if the value is in $x
 2618     my $contained = 2;
 2619     my $idx_contained = vsearch( $contained, $x, { mode => $mode } );
 2620     my $x_contained = $x->copy;
 2621     $x_contained->slice( $idx_contained ) .= 9;
 2622     
 2623     # if the value is not in $x
 2624     my $not_contained = 1.5;
 2625     my $idx_not_contained = vsearch( $not_contained, $x, { mode => $mode } );
 2626     my $x_not_contained = $x->copy;
 2627     $x_not_contained->slice( $idx_not_contained ) .= 9;
 2628     
 2629     print sprintf("%-23s%30s\n", '$x', $x);
 2630     print sprintf("%-23s%30s\n",   "$mode ($contained)", $x_contained);
 2631     print sprintf("%-23s%30s\n\n", "$mode ($not_contained)", $x_not_contained);
 2632   }
 2633   
 2634   # $x                     [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
 2635   # sample (2)             [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
 2636   # sample (1.5)           [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
 2637   # 
 2638   # $x                     [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
 2639   # insert_leftmost (2)    [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
 2640   # insert_leftmost (1.5)  [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
 2641   # 
 2642   # $x                     [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
 2643   # insert_rightmost (2)   [0 0 0 1 1 1 2 2 2 9 3 3 4 4 4]
 2644   # insert_rightmost (1.5) [0 0 0 1 1 1 9 2 2 3 3 3 4 4 4]
 2645   # 
 2646   # $x                     [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
 2647   # match (2)              [0 0 0 1 1 1 2 9 2 3 3 3 4 4 4]
 2648   # match (1.5)            [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
 2649   # 
 2650   # $x                     [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
 2651   # bin_inclusive (2)      [0 0 0 1 1 1 2 2 9 3 3 3 4 4 4]
 2652   # bin_inclusive (1.5)    [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
 2653   # 
 2654   # $x                     [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]
 2655   # bin_exclusive (2)      [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
 2656   # bin_exclusive (1.5)    [0 0 0 1 1 9 2 2 2 3 3 3 4 4 4]
 2657 
 2658 
 2659 Also see
 2660 L<B<vsearch_sample>|/vsearch_sample>,
 2661 L<B<vsearch_insert_leftmost>|/vsearch_insert_leftmost>,
 2662 L<B<vsearch_insert_rightmost>|/vsearch_insert_rightmost>,
 2663 L<B<vsearch_match>|/vsearch_match>,
 2664 L<B<vsearch_bin_inclusive>|/vsearch_bin_inclusive>, and
 2665 L<B<vsearch_bin_exclusive>|/vsearch_bin_exclusive>
 2666 
 2667 =cut
 2668 
 2669 sub vsearch {
 2670     my $opt = 'HASH' eq ref $_[-1]
 2671             ? pop
 2672         : { mode => 'sample' };
 2673 
 2674     croak( "unknown options to vsearch\n" )
 2675     if ( ! defined $opt->{mode} && keys %$opt )
 2676     || keys %$opt > 1;
 2677 
 2678     my $mode = $opt->{mode};
 2679     goto
 2680         $mode eq 'sample'           ? \&vsearch_sample
 2681       : $mode eq 'insert_leftmost'  ? \&vsearch_insert_leftmost
 2682       : $mode eq 'insert_rightmost' ? \&vsearch_insert_rightmost
 2683       : $mode eq 'match'            ? \&vsearch_match
 2684       : $mode eq 'bin_inclusive'    ? \&vsearch_bin_inclusive
 2685       : $mode eq 'bin_exclusive'    ? \&vsearch_bin_exclusive
 2686       :                               croak( "unknown vsearch mode: $mode\n" );
 2687 }
 2688 
 2689 *PDL::vsearch = \&vsearch;
 2690 #line 2691 "Primitive.pm"
 2691 
 2692 
 2693 
 2694 #line 948 "../../blib/lib/PDL/PP.pm"
 2695 
 2696 
 2697 
 2698 =head2 vsearch_sample
 2699 
 2700 =for sig
 2701 
 2702   Signature: (vals(); x(n); indx [o]idx())
 2703 
 2704 
 2705 =for ref
 2706 
 2707 Search for values in a sorted array, return index appropriate for sampling from a distribution
 2708 
 2709 =for usage
 2710 
 2711   $idx = vsearch_sample($vals, $x);
 2712 
 2713 C<$x> must be sorted, but may be in decreasing or increasing
 2714 order.
 2715 
 2716 
 2717 
 2718 B<vsearch_sample> returns an index I<I> for each value I<V> of C<$vals> appropriate
 2719 for sampling C<$vals>
 2720 
 2721 
 2722 
 2723                    
 2724 
 2725 I<I> has the following properties:
 2726 
 2727 =over
 2728 
 2729 =item *
 2730 
 2731 if C<$x> is sorted in increasing order
 2732 
 2733 
 2734           V <= x[0]  : I = 0
 2735   x[0]  < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
 2736   x[-1] < V          : I = $x->nelem -1
 2737 
 2738 
 2739 
 2740 =item *
 2741 
 2742 if C<$x> is sorted in decreasing order
 2743 
 2744 
 2745            V > x[0]  : I = 0
 2746   x[0]  >= V > x[-1] : I s.t. x[I] >= V > x[I+1]
 2747   x[-1] >= V         : I = $x->nelem - 1
 2748 
 2749 
 2750 
 2751 =back
 2752 
 2753 
 2754 
 2755 
 2756 If all elements of C<$x> are equal, I<< I = $x->nelem - 1 >>.
 2757 
 2758 If C<$x> contains duplicated elements, I<I> is the index of the
 2759 leftmost (by position in array) duplicate if I<V> matches.
 2760 
 2761 =for example
 2762 
 2763 This function is useful e.g. when you have a list of probabilities
 2764 for events and want to generate indices to events:
 2765 
 2766  $x = pdl(.01,.86,.93,1); # Barnsley IFS probabilities cumulatively
 2767  $y = random 20;
 2768  $c = vsearch_sample($y, $x); # Now, $c will have the appropriate distr.
 2769 
 2770 It is possible to use the L<cumusumover|PDL::Ufunc/cumusumover>
 2771 function to obtain cumulative probabilities from absolute probabilities.
 2772 
 2773 
 2774 
 2775 
 2776 
 2777 
 2778 
 2779 =for bad
 2780 
 2781 needs major (?) work to handles bad values
 2782 
 2783 =cut
 2784 #line 2785 "Primitive.pm"
 2785 
 2786 
 2787 
 2788 #line 950 "../../blib/lib/PDL/PP.pm"
 2789 
 2790 *vsearch_sample = \&PDL::vsearch_sample;
 2791 #line 2792 "Primitive.pm"
 2792 
 2793 
 2794 
 2795 #line 948 "../../blib/lib/PDL/PP.pm"
 2796 
 2797 
 2798 
 2799 =head2 vsearch_insert_leftmost
 2800 
 2801 =for sig
 2802 
 2803   Signature: (vals(); x(n); indx [o]idx())
 2804 
 2805 
 2806 =for ref
 2807 
 2808 Determine the insertion point for values in a sorted array, inserting before duplicates.
 2809 
 2810 =for usage
 2811 
 2812   $idx = vsearch_insert_leftmost($vals, $x);
 2813 
 2814 C<$x> must be sorted, but may be in decreasing or increasing
 2815 order.
 2816 
 2817 
 2818 
 2819 B<vsearch_insert_leftmost> returns an index I<I> for each value I<V> of
 2820 C<$vals> equal to the leftmost position (by index in array) within
 2821 C<$x> that I<V> may be inserted and still maintain the order in
 2822 C<$x>.
 2823 
 2824 Insertion at index I<I> involves shifting elements I<I> and higher of
 2825 C<$x> to the right by one and setting the now empty element at index
 2826 I<I> to I<V>.
 2827 
 2828 
 2829 
 2830                    
 2831 
 2832 I<I> has the following properties:
 2833 
 2834 =over
 2835 
 2836 =item *
 2837 
 2838 if C<$x> is sorted in increasing order
 2839 
 2840 
 2841           V <= x[0]  : I = 0
 2842   x[0]  < V <= x[-1] : I s.t. x[I-1] < V <= x[I]
 2843   x[-1] < V          : I = $x->nelem
 2844 
 2845 
 2846 
 2847 =item *
 2848 
 2849 if C<$x> is sorted in decreasing order
 2850 
 2851 
 2852            V >  x[0]  : I = -1
 2853   x[0]  >= V >= x[-1] : I s.t. x[I] >= V > x[I+1]
 2854   x[-1] >= V          : I = $x->nelem -1
 2855 
 2856 
 2857 
 2858 =back
 2859 
 2860 
 2861 
 2862 
 2863 If all elements of C<$x> are equal,
 2864 
 2865     i = 0
 2866 
 2867 If C<$x> contains duplicated elements, I<I> is the index of the
 2868 leftmost (by index in array) duplicate if I<V> matches.
 2869 
 2870 
 2871 
 2872 
 2873 
 2874 
 2875 
 2876 =for bad
 2877 
 2878 needs major (?) work to handles bad values
 2879 
 2880 =cut
 2881 #line 2882 "Primitive.pm"
 2882 
 2883 
 2884 
 2885 #line 950 "../../blib/lib/PDL/PP.pm"
 2886 
 2887 *vsearch_insert_leftmost = \&PDL::vsearch_insert_leftmost;
 2888 #line 2889 "Primitive.pm"
 2889 
 2890 
 2891 
 2892 #line 948 "../../blib/lib/PDL/PP.pm"
 2893 
 2894 
 2895 
 2896 =head2 vsearch_insert_rightmost
 2897 
 2898 =for sig
 2899 
 2900   Signature: (vals(); x(n); indx [o]idx())
 2901 
 2902 
 2903 =for ref
 2904 
 2905 Determine the insertion point for values in a sorted array, inserting after duplicates.
 2906 
 2907 =for usage
 2908 
 2909   $idx = vsearch_insert_rightmost($vals, $x);
 2910 
 2911 C<$x> must be sorted, but may be in decreasing or increasing
 2912 order.
 2913 
 2914 
 2915 
 2916 B<vsearch_insert_rightmost> returns an index I<I> for each value I<V> of
 2917 C<$vals> equal to the rightmost position (by index in array) within
 2918 C<$x> that I<V> may be inserted and still maintain the order in
 2919 C<$x>.
 2920 
 2921 Insertion at index I<I> involves shifting elements I<I> and higher of
 2922 C<$x> to the right by one and setting the now empty element at index
 2923 I<I> to I<V>.
 2924 
 2925 
 2926 
 2927                    
 2928 
 2929 I<I> has the following properties:
 2930 
 2931 =over
 2932 
 2933 =item *
 2934 
 2935 if C<$x> is sorted in increasing order
 2936 
 2937 
 2938            V < x[0]  : I = 0
 2939   x[0]  <= V < x[-1] : I s.t. x[I-1] <= V < x[I]
 2940   x[-1] <= V         : I = $x->nelem
 2941 
 2942 
 2943 
 2944 =item *
 2945 
 2946 if C<$x> is sorted in decreasing order
 2947 
 2948 
 2949           V >= x[0]  : I = -1
 2950   x[0]  > V >= x[-1] : I s.t. x[I] >= V > x[I+1]
 2951   x[-1] > V          : I = $x->nelem -1
 2952 
 2953 
 2954 
 2955 =back
 2956 
 2957 
 2958 
 2959 
 2960 If all elements of C<$x> are equal,
 2961 
 2962     i = $x->nelem - 1
 2963 
 2964 If C<$x> contains duplicated elements, I<I> is the index of the
 2965 leftmost (by index in array) duplicate if I<V> matches.
 2966 
 2967 
 2968 
 2969 
 2970 
 2971 
 2972 
 2973 =for bad
 2974 
 2975 needs major (?) work to handles bad values
 2976 
 2977 =cut
 2978 #line 2979 "Primitive.pm"
 2979 
 2980 
 2981 
 2982 #line 950 "../../blib/lib/PDL/PP.pm"
 2983 
 2984 *vsearch_insert_rightmost = \&PDL::vsearch_insert_rightmost;
 2985 #line 2986 "Primitive.pm"
 2986 
 2987 
 2988 
 2989 #line 948 "../../blib/lib/PDL/PP.pm"
 2990 
 2991 
 2992 
 2993 =head2 vsearch_match
 2994 
 2995 =for sig
 2996 
 2997   Signature: (vals(); x(n); indx [o]idx())
 2998 
 2999 
 3000 =for ref
 3001 
 3002 Match values against a sorted array.
 3003 
 3004 =for usage
 3005 
 3006   $idx = vsearch_match($vals, $x);
 3007 
 3008 C<$x> must be sorted, but may be in decreasing or increasing
 3009 order.
 3010 
 3011 
 3012 
 3013 B<vsearch_match> returns an index I<I> for each value I<V> of
 3014 C<$vals>.  If I<V> matches an element in C<$x>, I<I> is the
 3015 index of that element, otherwise it is I<-( insertion_point + 1 )>,
 3016 where I<insertion_point> is an index in C<$x> where I<V> may be
 3017 inserted while maintaining the order in C<$x>.  If C<$x> has
 3018 duplicated values, I<I> may refer to any of them.
 3019 
 3020 
 3021 
 3022                    
 3023 
 3024 
 3025 
 3026 
 3027 
 3028 =for bad
 3029 
 3030 needs major (?) work to handles bad values
 3031 
 3032 =cut
 3033 #line 3034 "Primitive.pm"
 3034 
 3035 
 3036 
 3037 #line 950 "../../blib/lib/PDL/PP.pm"
 3038 
 3039 *vsearch_match = \&PDL::vsearch_match;
 3040 #line 3041 "Primitive.pm"
 3041 
 3042 
 3043 
 3044 #line 948 "../../blib/lib/PDL/PP.pm"
 3045 
 3046 
 3047 
 3048 =head2 vsearch_bin_inclusive
 3049 
 3050 =for sig
 3051 
 3052   Signature: (vals(); x(n); indx [o]idx())
 3053 
 3054 
 3055 =for ref
 3056 
 3057 Determine the index for values in a sorted array of bins, lower bound inclusive.
 3058 
 3059 =for usage
 3060 
 3061   $idx = vsearch_bin_inclusive($vals, $x);
 3062 
 3063 C<$x> must be sorted, but may be in decreasing or increasing
 3064 order.
 3065 
 3066 
 3067 
 3068 C<$x> represents the edges of contiguous bins, with the first and
 3069 last elements representing the outer edges of the outer bins, and the
 3070 inner elements the shared bin edges.
 3071 
 3072 The lower bound of a bin is inclusive to the bin, its outer bound is exclusive to it.
 3073 B<vsearch_bin_inclusive> returns an index I<I> for each value I<V> of C<$vals>
 3074 
 3075 
 3076 
 3077                    
 3078 
 3079 I<I> has the following properties:
 3080 
 3081 =over
 3082 
 3083 =item *
 3084 
 3085 if C<$x> is sorted in increasing order
 3086 
 3087 
 3088            V < x[0]  : I = -1
 3089   x[0]  <= V < x[-1] : I s.t. x[I] <= V < x[I+1]
 3090   x[-1] <= V         : I = $x->nelem - 1
 3091 
 3092 
 3093 
 3094 =item *
 3095 
 3096 if C<$x> is sorted in decreasing order
 3097 
 3098 
 3099            V >= x[0]  : I = 0
 3100   x[0]  >  V >= x[-1] : I s.t. x[I+1] > V >= x[I]
 3101   x[-1] >  V          : I = $x->nelem
 3102 
 3103 
 3104 
 3105 =back
 3106 
 3107 
 3108 
 3109 
 3110 If all elements of C<$x> are equal,
 3111 
 3112     i = $x->nelem - 1
 3113 
 3114 If C<$x> contains duplicated elements, I<I> is the index of the
 3115 righmost (by index in array) duplicate if I<V> matches.
 3116 
 3117 
 3118 
 3119 
 3120 
 3121 
 3122 
 3123 =for bad
 3124 
 3125 needs major (?) work to handles bad values
 3126 
 3127 =cut
 3128 #line 3129 "Primitive.pm"
 3129 
 3130 
 3131 
 3132 #line 950 "../../blib/lib/PDL/PP.pm"
 3133 
 3134 *vsearch_bin_inclusive = \&PDL::vsearch_bin_inclusive;
 3135 #line 3136 "Primitive.pm"
 3136 
 3137 
 3138 
 3139 #line 948 "../../blib/lib/PDL/PP.pm"
 3140 
 3141 
 3142 
 3143 =head2 vsearch_bin_exclusive
 3144 
 3145 =for sig
 3146 
 3147   Signature: (vals(); x(n); indx [o]idx())
 3148 
 3149 
 3150 =for ref
 3151 
 3152 Determine the index for values in a sorted array of bins, lower bound exclusive.
 3153 
 3154 =for usage
 3155 
 3156   $idx = vsearch_bin_exclusive($vals, $x);
 3157 
 3158 C<$x> must be sorted, but may be in decreasing or increasing
 3159 order.
 3160 
 3161 
 3162 
 3163 C<$x> represents the edges of contiguous bins, with the first and
 3164 last elements representing the outer edges of the outer bins, and the
 3165 inner elements the shared bin edges.
 3166 
 3167 The lower bound of a bin is exclusive to the bin, its upper bound is inclusive to it.
 3168 B<vsearch_bin_exclusive> returns an index I<I> for each value I<V> of C<$vals>.
 3169 
 3170 
 3171 
 3172                    
 3173 
 3174 I<I> has the following properties:
 3175 
 3176 =over
 3177 
 3178 =item *
 3179 
 3180 if C<$x> is sorted in increasing order
 3181 
 3182 
 3183            V <= x[0]  : I = -1
 3184   x[0]  <  V <= x[-1] : I s.t. x[I] < V <= x[I+1]
 3185   x[-1] <  V          : I = $x->nelem - 1
 3186 
 3187 
 3188 
 3189 =item *
 3190 
 3191 if C<$x> is sorted in decreasing order
 3192 
 3193 
 3194            V >  x[0]  : I = 0
 3195   x[0]  >= V >  x[-1] : I s.t. x[I-1] >= V > x[I]
 3196   x[-1] >= V          : I = $x->nelem
 3197 
 3198 
 3199 
 3200 =back
 3201 
 3202 
 3203 
 3204 
 3205 If all elements of C<$x> are equal,
 3206 
 3207     i = $x->nelem - 1
 3208 
 3209 If C<$x> contains duplicated elements, I<I> is the index of the
 3210 righmost (by index in array) duplicate if I<V> matches.
 3211 
 3212 
 3213 
 3214 
 3215 
 3216 
 3217 
 3218 =for bad
 3219 
 3220 needs major (?) work to handles bad values
 3221 
 3222 =cut
 3223 #line 3224 "Primitive.pm"
 3224 
 3225 
 3226 
 3227 #line 950 "../../blib/lib/PDL/PP.pm"
 3228 
 3229 *vsearch_bin_exclusive = \&PDL::vsearch_bin_exclusive;
 3230 #line 3231 "Primitive.pm"
 3231 
 3232 
 3233 
 3234 #line 948 "../../blib/lib/PDL/PP.pm"
 3235 
 3236 
 3237 
 3238 =head2 interpolate
 3239 
 3240 =for sig
 3241 
 3242   Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
 3243 
 3244 
 3245 =for ref
 3246 
 3247 routine for 1D linear interpolation
 3248 
 3249 =for usage
 3250 
 3251  ( $yi, $err ) = interpolate($xi, $x, $y)
 3252 
 3253 Given a set of points C<($x,$y)>, use linear interpolation
 3254 to find the values C<$yi> at a set of points C<$xi>.
 3255 
 3256 C<interpolate> uses a binary search to find the suspects, er...,
 3257 interpolation indices and therefore abscissas (ie C<$x>)
 3258 have to be I<strictly> ordered (increasing or decreasing).
 3259 For interpolation at lots of
 3260 closely spaced abscissas an approach that uses the last index found as
 3261 a start for the next search can be faster (compare Numerical Recipes
 3262 C<hunt> routine). Feel free to implement that on top of the binary
 3263 search if you like. For out of bounds values it just does a linear
 3264 extrapolation and sets the corresponding element of C<$err> to 1,
 3265 which is otherwise 0.
 3266 
 3267 See also L</interpol>, which uses the same routine,
 3268 differing only in the handling of extrapolation - an error message
 3269 is printed rather than returning an error ndarray.
 3270 
 3271 
 3272 
 3273 =for bad
 3274 
 3275 needs major (?) work to handles bad values
 3276 
 3277 =cut
 3278 #line 3279 "Primitive.pm"
 3279 
 3280 
 3281 
 3282 #line 950 "../../blib/lib/PDL/PP.pm"
 3283 
 3284 *interpolate = \&PDL::interpolate;
 3285 #line 3286 "Primitive.pm"
 3286 
 3287 
 3288 
 3289 #line 3029 "primitive.pd"
 3290 
 3291 =head2 interpol
 3292 
 3293 =for sig
 3294 
 3295  Signature: (xi(); x(n); y(n); [o] yi())
 3296 
 3297 =for ref
 3298 
 3299 routine for 1D linear interpolation
 3300 
 3301 =for usage
 3302 
 3303  $yi = interpol($xi, $x, $y)
 3304 
 3305 C<interpol> uses the same search method as L</interpolate>,
 3306 hence C<$x> must be I<strictly> ordered (either increasing or decreasing).
 3307 The difference occurs in the handling of out-of-bounds values; here
 3308 an error message is printed.
 3309 
 3310 =cut
 3311 
 3312 # kept in for backwards compatability
 3313 sub interpol ($$$;$) {
 3314     my $xi = shift;
 3315     my $x  = shift;
 3316     my $y  = shift;
 3317     my $yi = @_ == 1 ? $_[0] : PDL->null;
 3318     interpolate( $xi, $x, $y, $yi, my $err = PDL->null );
 3319     print "some values had to be extrapolated\n"
 3320     if any $err;
 3321     return $yi if @_ == 0;
 3322 } # sub: interpol()
 3323 *PDL::interpol = \&interpol;
 3324 #line 3325 "Primitive.pm"
 3325 
 3326 
 3327 
 3328 #line 3067 "primitive.pd"
 3329 
 3330 =head2 interpND
 3331 
 3332 =for ref
 3333 
 3334 Interpolate values from an N-D ndarray, with switchable method
 3335 
 3336 =for example
 3337 
 3338   $source = 10*xvals(10,10) + yvals(10,10);
 3339   $index = pdl([[2.2,3.5],[4.1,5.0]],[[6.0,7.4],[8,9]]);
 3340   print $source->interpND( $index );
 3341 
 3342 InterpND acts like L<indexND|PDL::Slices/indexND>,
 3343 collapsing C<$index> by lookup
 3344 into C<$source>; but it does interpolation rather than direct sampling.
 3345 The interpolation method and boundary condition are switchable via
 3346 an options hash.
 3347 
 3348 By default, linear or sample interpolation is used, with constant
 3349 value outside the boundaries of the source pdl.  No dataflow occurs,
 3350 because in general the output is computed rather than indexed.
 3351 
 3352 All the interpolation methods treat the pixels as value-centered, so
 3353 the C<sample> method will return C<< $a->(0) >> for coordinate values on
 3354 the set [-0.5,0.5), and all methods will return C<< $a->(1) >> for
 3355 a coordinate value of exactly 1.
 3356 
 3357 
 3358 Recognized options:
 3359 
 3360 =over 3
 3361 
 3362 =item method
 3363 
 3364 Values can be:
 3365 
 3366 =over 3
 3367 
 3368 =item * 0, s, sample, Sample (default for integer source types)
 3369 
 3370 The nearest value is taken. Pixels are regarded as centered on their
 3371 respective integer coordinates (no offset from the linear case).
 3372 
 3373 =item * 1, l, linear, Linear (default for floating point source types)
 3374 
 3375 The values are N-linearly interpolated from an N-dimensional cube of size 2.
 3376 
 3377 =item * 3, c, cube, cubic, Cubic
 3378 
 3379 The values are interpolated using a local cubic fit to the data.  The
 3380 fit is constrained to match the original data and its derivative at the
 3381 data points.  The second derivative of the fit is not continuous at the
 3382 data points.  Multidimensional datasets are interpolated by the
 3383 successive-collapse method.
 3384 
 3385 (Note that the constraint on the first derivative causes a small amount
 3386 of ringing around sudden features such as step functions).
 3387 
 3388 =item * f, fft, fourier, Fourier
 3389 
 3390 The source is Fourier transformed, and the interpolated values are
 3391 explicitly calculated from the coefficients.  The boundary condition
 3392 option is ignored -- periodic boundaries are imposed.
 3393 
 3394 If you pass in the option "fft", and it is a list (ARRAY) ref, then it
 3395 is a stash for the magnitude and phase of the source FFT.  If the list
 3396 has two elements then they are taken as already computed; otherwise
 3397 they are calculated and put in the stash.
 3398 
 3399 =back
 3400 
 3401 =item b, bound, boundary, Boundary
 3402 
 3403 This option is passed unmodified into L<indexND|PDL::Slices/indexND>,
 3404 which is used as the indexing engine for the interpolation.
 3405 Some current allowed values are 'extend', 'periodic', 'truncate', and 'mirror'
 3406 (default is 'truncate').
 3407 
 3408 =item bad
 3409 
 3410 contains the fill value used for 'truncate' boundary.  (default 0)
 3411 
 3412 =item fft
 3413 
 3414 An array ref whose associated list is used to stash the FFT of the source
 3415 data, for the FFT method.
 3416 
 3417 =back
 3418 
 3419 =cut
 3420 
 3421 *interpND = *PDL::interpND;
 3422 sub PDL::interpND {
 3423   my $source = shift;
 3424   my $index = shift;
 3425   my $options = shift;
 3426 
 3427   barf 'Usage: interp_nd($source,$index,[{%options}])\n'
 3428     if(defined $options   and    ref $options ne 'HASH');
 3429 
 3430   my($opt) = (defined $options) ? $options : {};
 3431 
 3432   my($method)   = $opt->{m} || $opt->{meth} || $opt->{method} || $opt->{Method};
 3433   $method //= $source->type->integer ? 'sample' : 'linear';
 3434 
 3435   my($boundary) = $opt->{b} || $opt->{boundary} || $opt->{Boundary} || $opt->{bound} || $opt->{Bound} || 'extend';
 3436   my($bad) = $opt->{bad} || $opt->{Bad} || 0.0;
 3437 
 3438   if($method =~ m/^s(am(p(le)?)?)?/i) {
 3439     return $source->range(PDL::Math::floor($index+0.5),0,$boundary);
 3440   }
 3441 
 3442   elsif (($method eq 1) || $method =~ m/^l(in(ear)?)?/i) {
 3443     ## key: (ith = index broadcast; cth = cube broadcast; sth = source broadcast)
 3444     my $d = $index->dim(0);
 3445     my $di = $index->ndims - 1;
 3446 
 3447     # Grab a 2-on-a-side n-cube around each desired pixel
 3448     my $samp = $source->range($index->floor,2,$boundary); # (ith, cth, sth)
 3449 
 3450     # Reorder to put the cube dimensions in front and convert to a list
 3451     $samp = $samp->reorder( $di .. $di+$d-1,
 3452                 0 .. $di-1,
 3453                 $di+$d .. $samp->ndims-1) # (cth, ith, sth)
 3454                   ->clump($d); # (clst, ith, sth)
 3455 
 3456     # Enumerate the corners of an n-cube and convert to a list
 3457     # (the 'x' is the normal perl repeat operator)
 3458     my $crnr = PDL::Basic::ndcoords( (2) x $index->dim(0) ) # (index,cth)
 3459              ->mv(0,-1)->clump($index->dim(0))->mv(-1,0); # (index, clst)
 3460 
 3461     # a & b are the weighting coefficients.
 3462     my($x,$y);
 3463     my($indexwhere);
 3464     ($indexwhere = $index->where( 0 * $index )) .= -10; # Change NaN to invalid
 3465     {
 3466       my $bb = PDL::Math::floor($index);
 3467       $x = ($index - $bb)     -> dummy(1,$crnr->dim(1)); # index, clst, ith
 3468       $y = ($bb + 1 - $index) -> dummy(1,$crnr->dim(1)); # index, clst, ith
 3469     }
 3470 
 3471     # Use 1/0 corners to select which multiplier happens, multiply
 3472     # 'em all together to get sample weights, and sum to get the answer.
 3473     my $out0 =  ( ($x * ($crnr==1) + $y * ($crnr==0)) #index, clst, ith
 3474          -> prodover                          #clst, ith
 3475          );
 3476 
 3477     my $out = ($out0 * $samp)->sumover; # ith, sth
 3478 
 3479     # Work around BAD-not-being-contagious bug in PDL <= 2.6 bad handling code  --CED 3-April-2013
 3480     if ($source->badflag) {
 3481     my $baddies = $samp->isbad->orover;
 3482     $out = $out->setbadif($baddies);
 3483     }
 3484 
 3485     return $out;
 3486 
 3487   } elsif(($method eq 3) || $method =~ m/^c(u(b(e|ic)?)?)?/i) {
 3488 
 3489       my ($d,@di) = $index->dims;
 3490       my $di = $index->ndims - 1;
 3491 
 3492       # Grab a 4-on-a-side n-cube around each desired pixel
 3493       my $samp = $source->range($index->floor - 1,4,$boundary) #ith, cth, sth
 3494       ->reorder( $di .. $di+$d-1, 0..$di-1, $di+$d .. $source->ndims-1 );
 3495                        # (cth, ith, sth)
 3496 
 3497       # Make a cube of the subpixel offsets, and expand its dims to
 3498       # a 4-on-a-side N-1 cube, to match the slices of $samp (used below).
 3499       my $y = $index - $index->floor;
 3500       for my $i(1..$d-1) {
 3501       $y = $y->dummy($i,4);
 3502       }
 3503 
 3504       # Collapse by interpolation, one dimension at a time...
 3505       for my $i(0..$d-1) {
 3506       my $a0 = $samp->slice("(1)");    # Just-under-sample
 3507       my $a1 = $samp->slice("(2)");    # Just-over-sample
 3508       my $a1a0 = $a1 - $a0;
 3509 
 3510       my $gradient = 0.5 * ($samp->slice("2:3")-$samp->slice("0:1"));
 3511       my $s0 = $gradient->slice("(0)");   # Just-under-gradient
 3512       my $s1 = $gradient->slice("(1)");   # Just-over-gradient
 3513 
 3514       my $bb = $y->slice("($i)");
 3515 
 3516       # Collapse the sample...
 3517       $samp = ( $a0 +
 3518             $bb * (
 3519                $s0  +
 3520                $bb * ( (3 * $a1a0 - 2*$s0 - $s1) +
 3521                    $bb * ( $s1 + $s0 - 2*$a1a0 )
 3522                    )
 3523                )
 3524             );
 3525 
 3526       # "Collapse" the subpixel offset...
 3527       $y = $y->slice(":,($i)");
 3528       }
 3529 
 3530       return $samp;
 3531 
 3532   } elsif($method =~ m/^f(ft|ourier)?/i) {
 3533 
 3534      local $@;
 3535      eval "use PDL::FFT;";
 3536      my $fftref = $opt->{fft};
 3537      $fftref = [] unless(ref $fftref eq 'ARRAY');
 3538      if(@$fftref != 2) {
 3539      my $x = $source->copy;
 3540      my $y = zeroes($source);
 3541      fftnd($x,$y);
 3542      $fftref->[0] = sqrt($x*$x+$y*$y) / $x->nelem;
 3543      $fftref->[1] = - atan2($y,$x);
 3544      }
 3545 
 3546      my $i;
 3547      my $c = PDL::Basic::ndcoords($source);               # (dim, source-dims)
 3548      for $i(1..$index->ndims-1) {
 3549      $c = $c->dummy($i,$index->dim($i))
 3550      }
 3551      my $id = $index->ndims-1;
 3552      my $phase = (($c * $index * 3.14159 * 2 / pdl($source->dims))
 3553           ->sumover) # (index-dims, source-dims)
 3554               ->reorder($id..$id+$source->ndims-1,0..$id-1); # (src, index)
 3555 
 3556      my $phref = $fftref->[1]->copy;        # (source-dims)
 3557      my $mag = $fftref->[0]->copy;          # (source-dims)
 3558 
 3559      for $i(1..$index->ndims-1) {
 3560      $phref = $phref->dummy(-1,$index->dim($i));
 3561      $mag = $mag->dummy(-1,$index->dim($i));
 3562      }
 3563      my $out = cos($phase + $phref ) * $mag;
 3564      $out = $out->clump($source->ndims)->sumover;
 3565 
 3566      return $out;
 3567  }  else {
 3568      barf("interpND: unknown method '$method'; valid ones are 'linear' and 'sample'.\n");
 3569  }
 3570 }
 3571 #line 3572 "Primitive.pm"
 3572 
 3573 
 3574 
 3575 #line 3317 "primitive.pd"
 3576 
 3577 =head2 one2nd
 3578 
 3579 =for ref
 3580 
 3581 Converts a one dimensional index ndarray to a set of ND coordinates
 3582 
 3583 =for usage
 3584 
 3585  @coords=one2nd($x, $indices)
 3586 
 3587 returns an array of ndarrays containing the ND indexes corresponding to
 3588 the one dimensional list indices. The indices are assumed to
 3589 correspond to array C<$x> clumped using C<clump(-1)>. This routine is
 3590 used in the old vector form of L</whichND>, but is useful on
 3591 its own occasionally.
 3592 
 3593 Returned ndarrays have the L<indx|PDL::Core/indx> datatype.  C<$indices> can have
 3594 values larger than C<< $x->nelem >> but negative values in C<$indices>
 3595 will not give the answer you expect.
 3596 
 3597 =for example
 3598 
 3599  pdl> $x=pdl [[[1,2],[-1,1]], [[0,-3],[3,2]]]; $c=$x->clump(-1)
 3600  pdl> $maxind=maximum_ind($c); p $maxind;
 3601  6
 3602  pdl> print one2nd($x, maximum_ind($c))
 3603  0 1 1
 3604  pdl> p $x->at(0,1,1)
 3605  3
 3606 
 3607 =cut
 3608 
 3609 *one2nd = \&PDL::one2nd;
 3610 sub PDL::one2nd {
 3611   barf "Usage: one2nd \$array \$indices\n" if @_ != 2;
 3612   my ($x, $ind)=@_;
 3613   my @dimension=$x->dims;
 3614   $ind = indx($ind);
 3615   my(@index);
 3616   my $count=0;
 3617   foreach (@dimension) {
 3618     $index[$count++]=$ind % $_;
 3619     $ind /= $_;
 3620   }
 3621   return @index;
 3622 }
 3623 #line 3624 "Primitive.pm"
 3624 
 3625 
 3626 
 3627 #line 948 "../../blib/lib/PDL/PP.pm"
 3628 
 3629 
 3630 
 3631 =head2 which
 3632 
 3633 =for sig
 3634 
 3635   Signature: (mask(n); indx [o] inds(n); indx [o]lastout())
 3636 
 3637 
 3638 =for ref
 3639 
 3640 Returns indices of non-zero values from a 1-D PDL
 3641 
 3642 =for usage
 3643 
 3644  $i = which($mask);
 3645 
 3646 returns a pdl with indices for all those elements that are nonzero in
 3647 the mask. Note that the returned indices will be 1D. If you feed in a
 3648 multidimensional mask, it will be flattened before the indices are
 3649 calculated.  See also L</whichND> for multidimensional masks.
 3650 
 3651 If you want to index into the original mask or a similar ndarray
 3652 with output from C<which>, remember to flatten it before calling index:
 3653 
 3654   $data = random 5, 5;
 3655   $idx = which $data > 0.5; # $idx is now 1D
 3656   $bigsum = $data->flat->index($idx)->sum;  # flatten before indexing
 3657 
 3658 Compare also L</where> for similar functionality.
 3659 
 3660 SEE ALSO:
 3661 
 3662 L</which_both> returns separately the indices of both nonzero and zero
 3663 values in the mask.
 3664 
 3665 L</where_both> returns separately slices of both nonzero and zero
 3666 values in the mask.
 3667 
 3668 L</where> returns associated values from a data PDL, rather than
 3669 indices into the mask PDL.
 3670 
 3671 L</whichND> returns N-D indices into a multidimensional PDL.
 3672 
 3673 =for example
 3674 
 3675  pdl> $x = sequence(10); p $x
 3676  [0 1 2 3 4 5 6 7 8 9]
 3677  pdl> $indx = which($x>6); p $indx
 3678  [7 8 9]
 3679 
 3680 
 3681 
 3682 =for bad
 3683 
 3684 which processes bad values.
 3685 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 3686 
 3687 
 3688 =cut
 3689 #line 3690 "Primitive.pm"
 3690 
 3691 
 3692 
 3693 #line 949 "../../blib/lib/PDL/PP.pm"
 3694 
 3695    sub which { my ($this,$out) = @_;
 3696         $this = $this->flat;
 3697         $out //= $this->nullcreate;
 3698         PDL::_which_int($this,$out,my $lastout = $this->nullcreate);
 3699         my $lastoutmax = $lastout->max->sclr;
 3700         $lastoutmax ? $out->slice('0:'.($lastoutmax-1))->sever : empty(indx);
 3701    }
 3702    *PDL::which = \&which;
 3703 #line 3704 "Primitive.pm"
 3704 
 3705 
 3706 
 3707 #line 950 "../../blib/lib/PDL/PP.pm"
 3708 
 3709 *which = \&PDL::which;
 3710 #line 3711 "Primitive.pm"
 3711 
 3712 
 3713 
 3714 #line 948 "../../blib/lib/PDL/PP.pm"
 3715 
 3716 
 3717 
 3718 =head2 which_both
 3719 
 3720 =for sig
 3721 
 3722   Signature: (mask(n); indx [o] inds(n); indx [o]notinds(n); indx [o]lastout(); indx [o]lastoutn())
 3723 
 3724 
 3725 =for ref
 3726 
 3727 Returns indices of nonzero and zero values in a mask PDL
 3728 
 3729 =for usage
 3730 
 3731  ($i, $c_i) = which_both($mask);
 3732 
 3733 This works just as L</which>, but the complement of C<$i> will be in
 3734 C<$c_i>.
 3735 
 3736 =for example
 3737 
 3738  pdl> p $x = sequence(10)
 3739  [0 1 2 3 4 5 6 7 8 9]
 3740  pdl> ($big, $small) = which_both($x >= 5); p "$big\n$small"
 3741  [5 6 7 8 9]
 3742  [0 1 2 3 4]
 3743 
 3744 
 3745 
 3746 =for bad
 3747 
 3748 which_both processes bad values.
 3749 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
 3750 
 3751 
 3752 =cut
 3753 #line 3754 "Primitive.pm"
 3754 
 3755 
 3756 
 3757 #line 949 "../../blib/lib/PDL/PP.pm"
 3758 
 3759    sub which_both { my ($this,$outi,$outni) = @_;
 3760         $this = $this->flat;
 3761         $outi //= $this->nullcreate;
 3762         $outni //= $this->nullcreate;
 3763         PDL::_which_both_int($this,$outi,$outni,my $lastout = $this->nullcreate,my $lastoutn = $this->nullcreate);
 3764         my $lastoutmax = $lastout->max->sclr;
 3765         $outi = $lastoutmax ? $outi->slice('0:'.($lastoutmax-1))->sever : empty(indx);
 3766         return $outi if !wantarray;
 3767         my $lastoutnmax = $lastoutn->max->sclr;
 3768         ($outi, $lastoutnmax ? $outni->slice('0:'.($lastoutnmax-1))->sever : empty(indx));
 3769    }
 3770    *PDL::which_both = \&which_both;
 3771 #line 3772 "Primitive.pm"
 3772 
 3773 
 3774 
 3775 #line 950 "../../blib/lib/PDL/PP.pm"
 3776 
 3777 *which_both = \&PDL::which_both;
 3778 #line 3779 "Primitive.pm"
 3779 
 3780 
 3781 
 3782 #line 3498 "primitive.pd"
 3783 
 3784 =head2 where
 3785 
 3786 =for ref
 3787 
 3788 Use a mask to select values from one or more data PDLs
 3789 
 3790 C<where> accepts one or more data ndarrays and a mask ndarray.  It
 3791 returns a list of output ndarrays, corresponding to the input data
 3792 ndarrays.  Each output ndarray is a 1-dimensional list of values in its
 3793 corresponding data ndarray. The values are drawn from locations where
 3794 the mask is nonzero.
 3795 
 3796 The output PDLs are still connected to the original data PDLs, for the
 3797 purpose of dataflow.
 3798 
 3799 C<where> combines the functionality of L</which> and L<index|PDL::Slices/index>
 3800 into a single operation.
 3801 
 3802 BUGS:
 3803 
 3804 While C<where> works OK for most N-dimensional cases, it does not
 3805 broadcast properly over (for example) the (N+1)th dimension in data
 3806 that is compared to an N-dimensional mask.  Use C<whereND> for that.
 3807 
 3808 =for usage
 3809 
 3810  $i = $x->where($x+5 > 0); # $i contains those elements of $x
 3811                            # where mask ($x+5 > 0) is 1
 3812  $i .= -5;  # Set those elements (of $x) to -5. Together, these
 3813             # commands clamp $x to a maximum of -5.
 3814 
 3815 It is also possible to use the same mask for several ndarrays with
 3816 the same call:
 3817 
 3818  ($i,$j,$k) = where($x,$y,$z, $x+5>0);
 3819 
 3820 Note: C<$i> is always 1-D, even if C<$x> is E<gt>1-D.
 3821 
 3822 WARNING: The first argument
 3823 (the values) and the second argument (the mask) currently have to have
 3824 the exact same dimensions (or horrible things happen). You *cannot*
 3825 broadcast over a smaller mask, for example.
 3826 
 3827 =cut
 3828 
 3829 sub PDL::where {
 3830     barf "Usage: where( \$pdl1, ..., \$pdlN, \$mask )\n" if @_ == 1;
 3831     if(@_ == 2) {
 3832     my($data,$mask) = @_;
 3833     $data = $_[0]->clump(-1) if $_[0]->getndims>1;
 3834     $mask = $_[1]->clump(-1) if $_[0]->getndims>1;
 3835     return $data->index($mask->which());
 3836     } else {
 3837     if($_[-1]->getndims > 1) {
 3838         my $mask = $_[-1]->clump(-1)->which;
 3839         return map {$_->clump(-1)->index($mask)} @_[0..$#_-1];
 3840     } else {
 3841         my $mask = $_[-1]->which;
 3842         return map {$_->index($mask)} @_[0..$#_-1];
 3843     }
 3844     }
 3845 }
 3846 *where = \&PDL::where;
 3847 #line 3848 "Primitive.pm"
 3848 
 3849 
 3850 
 3851 #line 3566 "primitive.pd"
 3852 
 3853 =head2 where_both
 3854 
 3855 =for ref
 3856 
 3857 Returns slices (non-zero in mask, zero) of an ndarray according to a mask
 3858 
 3859 =for usage
 3860 
 3861  ($match_vals, $non_match_vals) = where_both($pdl, $mask);
 3862 
 3863 This works like L</which_both>, but (flattened) data-flowing slices
 3864 rather than index-sets are returned.
 3865 
 3866 =for example
 3867 
 3868  pdl> p $x = sequence(10) + 2
 3869  [2 3 4 5 6 7 8 9 10 11]
 3870  pdl> ($big, $small) = where_both($x, $x > 5); p "$big\n$small"
 3871  [6 7 8 9 10 11]
 3872  [2 3 4 5]
 3873  pdl> p $big += 2, $small -= 1
 3874  [8 9 10 11 12 13] [1 2 3 4]
 3875  pdl> p $x
 3876  [1 2 3 4 8 9 10 11 12 13]
 3877 
 3878 =cut
 3879 
 3880 sub PDL::where_both {
 3881   barf "Usage: where_both(\$pdl, \$mask)\n" if @_ != 2;
 3882   my ($arr, $mask) = @_;  # $mask has 0==false, 1==true
 3883   my $arr_flat = $arr->clump(-1);
 3884   map $arr_flat->index1d($_), PDL::which_both($mask);
 3885 }
 3886 *where_both = \&PDL::where_both;
 3887 #line 3888 "Primitive.pm"
 3888 
 3889 
 3890 
 3891 #line 3604 "primitive.pd"
 3892 
 3893 =head2 whereND
 3894 
 3895 =for ref
 3896 
 3897 C<where> with support for ND masks and broadcasting
 3898 
 3899 C<whereND> accepts one or more data ndarrays and a
 3900 mask ndarray.  It returns a list of output ndarrays,
 3901 corresponding to the input data ndarrays.  The values
 3902 are drawn from locations where the mask is nonzero.
 3903 
 3904 C<whereND> differs from C<where> in that the mask
 3905 dimensionality is preserved which allows for
 3906 proper broadcasting of the selection operation over
 3907 higher dimensions.
 3908 
 3909 As with C<where> the output PDLs are still connected
 3910 to the original data PDLs, for the purpose of dataflow.
 3911 
 3912 =for usage
 3913 
 3914   $sdata = whereND $data, $mask
 3915   ($s1, $s2, ..., $sn) = whereND $d1, $d2, ..., $dn, $mask
 3916 
 3917   where
 3918 
 3919     $data is M dimensional
 3920     $mask is N < M dimensional
 3921     dims($data) 1..N == dims($mask) 1..N
 3922     with broadcasting over N+1 to M dimensions
 3923 
 3924 =for example
 3925 
 3926   $data   = sequence(4,3,2);   # example data array
 3927   $mask4  = (random(4)>0.5);   # example 1-D mask array, has $n4 true values
 3928   $mask43 = (random(4,3)>0.5); # example 2-D mask array, has $n43 true values
 3929   $sdat4  = whereND $data, $mask4;   # $sdat4 is a [$n4,3,2] pdl
 3930   $sdat43 = whereND $data, $mask43;  # $sdat43 is a [$n43,2] pdl
 3931 
 3932 Just as with C<where>, you can use the returned value in an
 3933 assignment. That means that both of these examples are valid:
 3934 
 3935   # Used to create a new slice stored in $sdat4:
 3936   $sdat4 = $data->whereND($mask4);
 3937   $sdat4 .= 0;
 3938   # Used in lvalue context:
 3939   $data->whereND($mask4) .= 0;
 3940 
 3941 SEE ALSO:
 3942 
 3943 L</whichND> returns N-D indices into a multidimensional PDL, from a mask.
 3944 
 3945 =cut
 3946 
 3947 sub PDL::whereND :lvalue {
 3948    barf "Usage: whereND( \$pdl1, ..., \$pdlN, \$mask )\n" if @_ == 1;
 3949    my $mask = pop @_;  # $mask has 0==false, 1==true
 3950    my @to_return;
 3951    my $n = PDL::sum($mask);
 3952    foreach my $arr (@_) {
 3953       my $sub_i = $mask * ones($arr);
 3954       my $where_sub_i = PDL::where($arr, $sub_i);
 3955       # count the number of dims in $mask and $arr
 3956       # $mask = a b c d e f.....
 3957       my @idims = dims($arr);
 3958       # ...and pop off the number of dims in $mask
 3959       foreach ( dims($mask) ) { shift(@idims) };
 3960       my $ndim = 0;
 3961       foreach my $id ($n, @idims[0..($#idims-1)]) {
 3962          $where_sub_i = $where_sub_i->splitdim($ndim++,$id) if $n>0;
 3963       }
 3964       push @to_return, $where_sub_i;
 3965    }
 3966    return (@to_return == 1) ? $to_return[0] : @to_return;
 3967 }
 3968 *whereND = \&PDL::whereND;
 3969 #line 3970 "Primitive.pm"
 3970 
 3971 
 3972 
 3973 #line 3685 "primitive.pd"
 3974 
 3975 =head2 whichND
 3976 
 3977 =for ref
 3978 
 3979 Return the coordinates of non-zero values in a mask.
 3980 
 3981 =for usage
 3982 
 3983 WhichND returns the N-dimensional coordinates of each nonzero value in
 3984 a mask PDL with any number of dimensions.  The returned values arrive
 3985 as an array-of-vectors suitable for use in
 3986 L<indexND|PDL::Slices/indexND> or L<range|PDL::Slices/range>.
 3987 
 3988  $coords = whichND($mask);
 3989 
 3990 returns a PDL containing the coordinates of the elements that are non-zero
 3991 in C<$mask>, suitable for use in L<PDL::Slices/indexND>. The 0th dimension contains the
 3992 full coordinate listing of each point; the 1st dimension lists all the points.
 3993 For example, if $mask has rank 4 and 100 matching elements, then $coords has
 3994 dimension 4x100.
 3995 
 3996 If no such elements exist, then whichND returns a structured empty PDL:
 3997 an Nx0 PDL that contains no values (but matches, broadcasting-wise, with
 3998 the vectors that would be produced if such elements existed).
 3999 
 4000 DEPRECATED BEHAVIOR IN LIST CONTEXT:
 4001 
 4002 whichND once delivered different values in list context than in scalar
 4003 context, for historical reasons.  In list context, it returned the
 4004 coordinates transposed, as a collection of 1-PDLs (one per dimension)
 4005 in a list.  This usage is deprecated in PDL 2.4.10, and will cause a
 4006 warning to be issued every time it is encountered.  To avoid the
 4007 warning, you can set the global variable "$PDL::whichND" to 's' to
 4008 get scalar behavior in all contexts, or to 'l' to get list behavior in
 4009 list context.
 4010 
 4011 In later versions of PDL, the deprecated behavior will disappear.  Deprecated
 4012 list context whichND expressions can be replaced with:
 4013 
 4014     @list = $x->whichND->mv(0,-1)->dog;
 4015 
 4016 SEE ALSO:
 4017 
 4018 L</which> finds coordinates of nonzero values in a 1-D mask.
 4019 
 4020 L</where> extracts values from a data PDL that are associated
 4021 with nonzero values in a mask PDL.
 4022 
 4023 L<PDL::Slices/indexND> can be fed the coordinates to return the values.
 4024 
 4025 =for example
 4026 
 4027  pdl> $s=sequence(10,10,3,4)
 4028  pdl> ($x, $y, $z, $w)=whichND($s == 203); p $x, $y, $z, $w
 4029  [3] [0] [2] [0]
 4030  pdl> print $s->at(list(cat($x,$y,$z,$w)))
 4031  203
 4032 
 4033 =cut
 4034 
 4035 *whichND = \&PDL::whichND;
 4036 sub PDL::whichND {
 4037   my $mask = PDL->topdl(shift);
 4038 
 4039   # List context: generate a perl list by dimension
 4040   if(wantarray) {
 4041       if(!defined($PDL::whichND)) {
 4042       printf STDERR "whichND: WARNING - list context deprecated. Set \$PDL::whichND. Details in pod.";
 4043       } elsif($PDL::whichND =~ m/l/i) {
 4044       # old list context enabled by setting $PDL::whichND to 'l'
 4045       my $ind=($mask->clump(-1))->which;
 4046       return $mask->one2nd($ind);
 4047       }
 4048       # if $PDL::whichND does not contain 'l' or 'L', fall through to scalar context
 4049   }
 4050 
 4051   # Scalar context: generate an N-D index ndarray
 4052   return PDL->new_from_specification(indx,$mask->ndims,0) if !$mask->nelem;
 4053   return $mask ? pdl(indx,0) : PDL->new_from_specification(indx,0)
 4054     if !$mask->getndims;
 4055 
 4056   my $ind = $mask->flat->which->dummy(0,$mask->getndims)->make_physical;
 4057   # In the empty case, explicitly return the correct type of structured empty
 4058   return PDL->new_from_specification(indx,$mask->ndims, 0) if !$ind->nelem;
 4059 
 4060   my $mult = ones(indx, $mask->getndims);
 4061   my @mdims = $mask->dims;
 4062 
 4063   for my $i (0..$#mdims-1) {
 4064    # use $tmp for 5.005_03 compatibility
 4065    (my $tmp = $mult->index($i+1)) .= $mult->index($i)*$mdims[$i];
 4066   }
 4067 
 4068   for my $i (0..$#mdims) {
 4069    my($s) = $ind->index($i);
 4070    $s /= $mult->index($i);
 4071    $s %= $mdims[$i];
 4072   }
 4073 
 4074   return $ind;
 4075 }
 4076 #line 4077 "Primitive.pm"
 4077 
 4078 
 4079 
 4080 #line 3794 "primitive.pd"
 4081 
 4082 =head2 setops
 4083 
 4084 =for ref
 4085 
 4086 Implements simple set operations like union and intersection
 4087 
 4088 =for usage
 4089 
 4090    Usage: $set = setops($x, <OPERATOR>, $y);
 4091 
 4092 The operator can be C<OR>, C<XOR> or C<AND>. This is then applied
 4093 to C<$x> viewed as a set and C<$y> viewed as a set. Set theory says
 4094 that a set may not have two or more identical elements, but setops
 4095 takes care of this for you, so C<$x=pdl(1,1,2)> is OK. The functioning
 4096 is as follows:
 4097 
 4098 =over
 4099 
 4100 =item C<OR>
 4101 
 4102 The resulting vector will contain the elements that are either in C<$x>
 4103 I<or> in C<$y> or both. This is the union in set operation terms
 4104 
 4105 =item C<XOR>
 4106 
 4107 The resulting vector will contain the elements that are either in C<$x>
 4108 or C<$y>, but not in both. This is
 4109 
 4110      Union($x, $y) - Intersection($x, $y)
 4111 
 4112 in set operation terms.
 4113 
 4114 =item C<AND>
 4115 
 4116 The resulting vector will contain the intersection of C<$x> and C<$y>, so
 4117 the elements that are in both C<$x> and C<$y>. Note that for convenience
 4118 this operation is also aliased to L</intersect>.
 4119 
 4120 =back
 4121 
 4122 It should be emphasized that these routines are used when one or both of
 4123 the sets C<$x>, C<$y> are hard to calculate or that you get from a separate
 4124 subroutine.
 4125 
 4126 Finally IDL users might be familiar with Craig Markwardt's C<cmset_op.pro>
 4127 routine which has inspired this routine although it was written independently
 4128 However the present routine has a few less options (but see the examples)
 4129 
 4130 =for example
 4131 
 4132 You will very often use these functions on an index vector, so that is
 4133 what we will show here. We will in fact something slightly silly. First
 4134 we will find all squares that are also cubes below 10000.
 4135 
 4136 Create a sequence vector:
 4137 
 4138   pdl> $x = sequence(10000)
 4139 
 4140 Find all odd and even elements:
 4141 
 4142   pdl> ($even, $odd) = which_both( ($x % 2) == 0)
 4143 
 4144 Find all squares
 4145 
 4146   pdl> $squares= which(ceil(sqrt($x)) == floor(sqrt($x)))
 4147 
 4148 Find all cubes (being careful with roundoff error!)
 4149 
 4150   pdl> $cubes= which(ceil($x**(1.0/3.0)) == floor($x**(1.0/3.0)+1e-6))
 4151 
 4152 Then find all squares that are cubes:
 4153 
 4154   pdl> $both = setops($squares, 'AND', $cubes)
 4155 
 4156 And print these (assumes that C<PDL::NiceSlice> is loaded!)
 4157 
 4158   pdl> p $x($both)
 4159    [0 1 64 729 4096]
 4160 
 4161 Then find all numbers that are either cubes or squares, but not both:
 4162 
 4163   pdl> $cube_xor_square = setops($squares, 'XOR', $cubes)
 4164 
 4165   pdl> p $cube_xor_square->nelem()
 4166    112
 4167 
 4168 So there are a total of 112 of these!
 4169 
 4170 Finally find all odd squares:
 4171 
 4172   pdl> $odd_squares = setops($squares, 'AND', $odd)
 4173 
 4174 
 4175 Another common occurrence is to want to get all objects that are
 4176 in C<$x> and in the complement of C<$y>. But it is almost always best
 4177 to create the complement explicitly since the universe that both are
 4178 taken from is not known. Thus use L</which_both> if possible
 4179 to keep track of complements.
 4180 
 4181 If this is impossible the best approach is to make a temporary:
 4182 
 4183 This creates an index vector the size of the universe of the sets and
 4184 set all elements in C<$y> to 0
 4185 
 4186   pdl> $tmp = ones($n_universe); $tmp($y) .= 0;
 4187 
 4188 This then finds the complement of C<$y>
 4189 
 4190   pdl> $C_b = which($tmp == 1);
 4191 
 4192 and this does the final selection:
 4193 
 4194   pdl> $set = setops($x, 'AND', $C_b)
 4195 
 4196 =cut
 4197 
 4198 *setops = \&PDL::setops;
 4199 
 4200 sub PDL::setops {
 4201 
 4202   my ($x, $op, $y)=@_;
 4203 
 4204   # Check that $x and $y are 1D.
 4205   if ($x->ndims() > 1 || $y->ndims() > 1) {
 4206      warn 'setops: $x and $y must be 1D - flattening them!'."\n";
 4207      $x = $x->flat;
 4208      $y = $y->flat;
 4209   }
 4210 
 4211   #Make sure there are no duplicate elements.
 4212   $x=$x->uniq;
 4213   $y=$y->uniq;
 4214 
 4215   my $result;
 4216 
 4217   if ($op eq 'OR') {
 4218     # Easy...
 4219     $result = uniq(append($x, $y));
 4220   } elsif ($op eq 'XOR') {
 4221     # Make ordered list of set union.
 4222     my $union = append($x, $y)->qsort;
 4223     # Index lists.
 4224     my $s1=zeroes(byte, $union->nelem());
 4225     my $s2=zeroes(byte, $union->nelem());
 4226 
 4227     # Find indices which are duplicated - these are to be excluded
 4228     #
 4229     # We do this by comparing x with x shifted each way.
 4230     my $i1 = which($union != rotate($union, 1));
 4231     my $i2 = which($union != rotate($union, -1));
 4232     #
 4233     # We then mark/mask these in the s1 and s2 arrays to indicate which ones
 4234     # are not equal to their neighbours.
 4235     #
 4236     my $ts;
 4237     ($ts = $s1->index($i1)) .= 1 if $i1->nelem() > 0;
 4238     ($ts = $s2->index($i2)) .= 1 if $i2->nelem() > 0;
 4239 
 4240     my $inds=which($s1 == $s2);
 4241 
 4242     if ($inds->nelem() > 0) {
 4243       return $union->index($inds);
 4244     } else {
 4245       return $inds;
 4246     }
 4247 
 4248   } elsif ($op eq 'AND') {
 4249     # The intersection of the arrays.
 4250 
 4251     # Make ordered list of set union.
 4252     my $union = append($x, $y)->qsort;
 4253 
 4254     return $union->where($union == rotate($union, -1));
 4255   } else {
 4256     print "The operation $op is not known!";
 4257     return -1;
 4258   }
 4259 
 4260 }
 4261 #line 4262 "Primitive.pm"
 4262 
 4263 
 4264 
 4265 #line 3977 "primitive.pd"
 4266 
 4267 =head2 intersect
 4268 
 4269 =for ref
 4270 
 4271 Calculate the intersection of two ndarrays
 4272 
 4273 =for usage
 4274 
 4275    Usage: $set = intersect($x, $y);
 4276 
 4277 This routine is merely a simple interface to L</setops>. See
 4278 that for more information
 4279 
 4280 =for example
 4281 
 4282 Find all numbers less that 100 that are of the form 2*y and 3*x
 4283 
 4284  pdl> $x=sequence(100)
 4285  pdl> $factor2 = which( ($x % 2) == 0)
 4286  pdl> $factor3 = which( ($x % 3) == 0)
 4287  pdl> $ii=intersect($factor2, $factor3)
 4288  pdl> p $x($ii)
 4289  [0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96]
 4290 
 4291 =cut
 4292 
 4293 *intersect = \&PDL::intersect;
 4294 
 4295 sub PDL::intersect {
 4296 
 4297    return setops($_[0], 'AND', $_[1]);
 4298 
 4299 }
 4300 #line 4301 "Primitive.pm"
 4301 
 4302 
 4303 
 4304 
 4305 
 4306 #line 4013 "primitive.pd"
 4307 
 4308 
 4309 =head1 AUTHOR
 4310 
 4311 Copyright (C) Tuomas J. Lukka 1997 (lukka@husc.harvard.edu). Contributions
 4312 by Christian Soeller (c.soeller@auckland.ac.nz), Karl Glazebrook
 4313 (kgb@aaoepp.aao.gov.au), Craig DeForest (deforest@boulder.swri.edu)
 4314 and Jarle Brinchmann (jarle@astro.up.pt)
 4315 All rights reserved. There is no warranty. You are allowed
 4316 to redistribute this software / documentation under certain
 4317 conditions. For details, see the file COPYING in the PDL
 4318 distribution. If this file is separated from the PDL distribution,
 4319 the copyright notice should be included in the file.
 4320 
 4321 Updated for CPAN viewing compatibility by David Mertens.
 4322 
 4323 =cut
 4324 #line 4325 "Primitive.pm"
 4325 
 4326 
 4327 
 4328 
 4329 # Exit with OK status
 4330 
 4331 1;