"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/t/primitive.t" (19 May 2022, 42533 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested text file into HTML format (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the latest Fossies "Diffs" side-by-side code changes report for "primitive.t": 2.079_vs_2.080.

    1 use strict;
    2 use warnings;
    3 use Test::More;
    4 use Test::Exception;
    5 use PDL::LiteF;
    6 use PDL::Types;
    7 
    8 sub tapprox {
    9     my($x,$y) = @_;
   10     $_ = pdl($_) for $x, $y;
   11     if(join(',',$x->dims) ne join(',',$y->dims)) {
   12        diag "APPROX: $x $y\n";
   13        diag "UNEQDIM\n";
   14        return 0;
   15     }
   16     my $d = max( abs($x-$y) );
   17     if($d >= 0.01) {
   18        diag "got=$x expected=$y\n";
   19     }
   20     $d < 0.01;
   21 }
   22 
   23 ok tapprox(pdl(1,2,3)->cmpvec(pdl(3,2,1)), -1), 'cmpvec less';
   24 ok tapprox(pdl(3,2,1)->cmpvec(pdl(1,2,3)), 1), 'cmpvec more';
   25 ok tapprox(pdl(3,2,1)->cmpvec(pdl(3,2,1)), 0), 'cmpvec same';
   26 is pdl('[1 BAD]')->cmpvec(pdl(3,2)).'', '-1', 'cmpvec bad before';
   27 is pdl('[BAD 1]')->cmpvec(pdl(3,2)).'', 'BAD', 'cmpvec bad';
   28 
   29 ok tapprox(pdl(3,2,1)->eqvec(pdl(1,2,3)), 0), 'eqvec diff';
   30 ok tapprox(pdl(3,2,1)->eqvec(pdl(3,2,1)), 1), 'eqvec same';
   31 is pdl('[2 1 BAD]')->eqvec(pdl(1,3,2)).'', 'BAD', 'eqvec bad before';
   32 is pdl('[2 BAD 1]')->eqvec(pdl(2,3,2)).'', 'BAD', 'eqvec bad';
   33 
   34 my $x = PDL->pdl([[5,4,3],[2,3,1.5]]);
   35 ok(tapprox($x->average(), PDL->pdl([4, 2.16666])), "average");
   36 ok(tapprox($x->sumover(), PDL->pdl([12, 6.5])), "sumover");
   37 ok(tapprox($x->prodover(), PDL->pdl([60, 9])), "prodover");
   38 
   39 my $y = PDL->pdl(4,3,1,0,0,0,0,5,2,0,3,6);
   40 my $c = ($y->xvals) + 10;
   41 ok(tapprox($y->where($y>4), PDL->pdl(5,6)), "where with >");
   42 ok(tapprox($y->which, PDL->pdl(0,1,2,7,8,10,11)), "which");
   43 ok(tapprox($c->where($y), PDL->pdl(10,11,12,17,18,20,21)), "where with mask");
   44 ok zeroes(3)->which->isempty, 'which on all-zero returns empty';
   45 
   46 $y = sequence(10) + 2;
   47 my ($big, $small) = where_both($y, $y > 5);
   48 $big += 2, $small -= 1;
   49 ok tapprox($big, pdl('[8 9 10 11 12 13]')), 'where_both big + 2 is right';
   50 ok tapprox($small, pdl('[1 2 3 4]')), 'where_both small - 2 is right';
   51 ok tapprox($y, pdl('[1 2 3 4 8 9 10 11 12 13]')), 'where_both dataflow affected orig';
   52 
   53 {
   54   my $orig = ones(byte, 300);
   55   my $xvals = $orig->xvals;
   56   is $xvals->at(280), 280, 'non-wrapped xvals from byte ndarray';
   57 }
   58 
   59 ##############################
   60 # originally in pptest
   61 $x = ones(byte,3000);
   62 dsumover($x,($y=null));
   63 is($y->get_datatype, $PDL_D, "get_datatype" );
   64 is($y->at, 3000, "at" );
   65 
   66 my $p = pdl [ 1, 2, 3, 4, 7, 9, 1, 1, 6, 2, 5];
   67 my $q = zeroes 5;
   68 minimum_n_ind $p, $q;
   69 ok(tapprox($q, pdl(0, 6, 7, 1, 9)), "minimum_n_ind");
   70 $q = minimum_n_ind($p, 5);
   71 ok(tapprox($q, pdl(0, 6, 7, 1, 9)), "minimum_n_ind usage 2");
   72 minimum_n_ind($p, $q = null, 5);
   73 ok(tapprox($q, pdl(0, 6, 7, 1, 9)), "minimum_n_ind usage 3");
   74 $p = pdl '[1 BAD 3 4 7 9 1 1 6 2 5]';
   75 $q = zeroes 5;
   76 minimum_n_ind $p, $q;
   77 is $q.'', '[0 6 7 9 2]', "minimum_n_ind BAD";
   78 $p = pdl '[1 BAD 3 4 BAD BAD]';
   79 $q = zeroes 5;
   80 minimum_n_ind $p, $q;
   81 is $q.'', '[0 2 3 BAD BAD]', "minimum_n_ind insufficient good";
   82 $p = pdl '[1 BAD 3 4 BAD BAD 3 1 5 8 9]';
   83 $q = zeroes 5;
   84 minimum_n_ind $p, $q;
   85 is $q.'', '[0 7 2 6 3]', "minimum_n_ind some bad, sufficient good";
   86 
   87 ##############################
   88 # check that our random functions work with Perl's srand
   89 TODO: {
   90    local $TODO = 'Some CPAN Testers fails for OpenBSD';
   91 
   92    srand 5;
   93    my $r1 = random 10;
   94    srand 5;
   95    my $r2 = random 10;
   96    ok(tapprox($r1, $r2), "random and srand");
   97 
   98    srand 10;
   99    $r1 = grandom 10;
  100    srand 10;
  101    $r2 = grandom 10;
  102    ok(tapprox($r1, $r2), "grandom and srand");
  103 }
  104 ##############################
  105 # Test that whichND works OK
  106 my $r = xvals(10,10)+10*yvals(10,10);
  107 $x = whichND( $r % 12 == 0 );
  108 
  109 # Nontrivial case gives correct coordinates
  110 ok(eval { sum($x != pdl([0,0],[2,1],[4,2],[6,3],[8,4],[0,6],[2,7],[4,8],[6,9]))==0 }, "whichND");
  111 is $x->type, 'indx', "whichND returns indx-type ndarray for non-trivial case";
  112 # Empty case gives matching Empty
  113 $x = whichND( $r*0 );
  114 is $x->nelem, 0, "whichND( 0*\$r ) gives an Empty PDL";
  115 is_deeply [$x->dims], [2,0], "whichND( 0*\$r ) is 2x0";
  116 is $x->type, 'indx', "whichND( 0*\$r) type is indx";
  117 
  118 # Scalar PDLs are treated as 1-PDLs
  119 $x = whichND(pdl(5));
  120 is $x->nelem, 1, "whichND scalar PDL";
  121 is $x, 0, "whichND scalar PDL";
  122 is $x->type, 'indx', "whichND returns indx ndarray for scalar ndarray mask";
  123 
  124 # Scalar empty case returns a 1-D vector of size 0
  125 $x = whichND(pdl(0));
  126 is $x->nelem, 0,  "whichND of 0 scalar is empty";
  127 is_deeply [$x->dims], [0], "whichND of 0 scalar: return 0 dim size is 0";
  128 is $x->type, 'indx', "whichND returns indx-type ndarray for scalar empty case";
  129 
  130 # Empty case returns Empty
  131 $y = whichND( which(pdl(0)) );                              
  132 is $y->nelem, 0, "whichND of Empty mask";
  133 is $y->type, 'indx', "whichND returns indx-type ndarray for empty case";
  134 
  135 # Nontrivial empty mask case returns matching Empty -- whichND(Empty[2x0x2]) should return Empty[3x0]
  136 $y = whichND(zeroes(2,0,2));
  137 is_deeply [$y->dims], [3,0], "whichND(Empty[2x0x2]) returns Empty[3x0]";
  138 
  139 $r = zeroes(7, 7); $r->set(3, 4, 1);
  140 is $r->whichND.'', <<EOF, 'whichND works right (was failing on 32-bit)';
  141 \n[\n [3 4]\n]
  142 EOF
  143 
  144 $x = pdl('[[i 2+3i] [4+5i 6+7i]]');
  145 ok all(approx $x->norm, pdl(<<'EOF')), 'native complex norm works' or diag $x->norm;
  146 [
  147  [0.267261i 0.534522+0.801783i]
  148  [0.356348+0.445435i 0.534522+0.623609i]
  149 ]
  150 EOF
  151 
  152 ##############################
  153 # Simple test case for interpND
  154 $x = xvals(10,10)+yvals(10,10)*10;
  155 my $index = cat(3+xvals(5,5)*0.25,7+yvals(5,5)*0.25)->reorder(2,0,1);
  156 my $z = 73+xvals(5,5)*0.25+2.5*yvals(5,5);
  157 eval { $y = $x->interpND($index) };
  158 is $@, '';
  159 is sum($y != $z), 0, "interpND";
  160 
  161 ##############################
  162 # Test glue
  163 $x = xvals(2,2,2);
  164 $y = yvals(2,2,2);
  165 $c = zvals(2,2,2);
  166 our $d;
  167 eval { $d = $x->glue(1,$y,$c) };
  168 is $@, '';
  169 ok(zcheck($d - pdl([[0,1],[0,1],[0,0],[1,1],[0,0],[0,0]],
  170                    [[0,1],[0,1],[0,0],[1,1],[1,1],[1,1]])), "glue");
  171 
  172 ##############################
  173 # test new empty ndarray handling
  174 $x = which ones(4) > 2;
  175 $y = $x->long;
  176 $c = $x->double;
  177 
  178 ok(isempty $x, "isempty");
  179 ok($y->avg == 0, "avg of Empty");
  180 ok(! any isfinite $c->average, "isfinite of Empty");
  181 
  182 ##############################
  183 # Test uniqvec
  184 $x = pdl([[0,1],[2,2],[0,1]]);
  185 $y = $x->uniqvec;
  186 eval { $c = all($y==pdl([[0,1],[2,2]])) };
  187 is $@, '';
  188 ok $c, "uniqvec";
  189 is $y->ndims, 2, "uniqvec";
  190 
  191 $x = pdl([[0,1]])->uniqvec;
  192 eval { $c = all($x==pdl([[0,1]])) };
  193 is $@, '';
  194 ok $c, "uniqvec";
  195 is $x->ndims, 2, "uniqvec";
  196 
  197 $x = pdl([[0,1,2]]); $x = $x->glue(1,$x,$x);
  198 $y = $x->uniqvec;
  199 eval { $c = all($y==pdl([0,1,2])) };
  200 is $@, '';
  201 ok $c, "uniqvec";
  202 is $y->ndims, 2, "uniqvec";
  203 
  204 ##############################
  205 # Test bad handling in selector
  206 $y = xvals(3);
  207 ok(tapprox($y->which,PDL->pdl(1,2)), "which");
  208 setbadat $y, 1;
  209 ok(tapprox($y->which,PDL->pdl([2])), "which w BAD");
  210 setbadat $y, 0;
  211 setbadat $y, 2;
  212 is($y->which->nelem,0, "which nelem w BAD");
  213 
  214 ############################
  215 # Test intersect & setops
  216 my $temp = sequence(10);
  217 $x = which(($temp % 2) == 0);
  218 $y = which(($temp % 3) == 0);
  219 $c = setops($x, 'AND', $y);
  220 ok(tapprox($c, pdl([0, 6])), "setops AND");
  221 ok(tapprox(intersect($x,$y),pdl([0,6])), "intersect same as setops AND");
  222 $c = setops($x,'OR',$y);
  223 ok(tapprox($c, pdl([0,2,3,4,6,8,9])), "setops OR");
  224 $c = setops($x,'XOR',$y);
  225 ok(tapprox($c, pdl([2,3,4,8,9])), "setops XOR");
  226 #Test intersect again
  227 my $intersect_test=intersect(pdl(1,-5,4,0), pdl(0,3,-5,2));
  228 ok (all($intersect_test==pdl(-5,0)), 'Intersect test values');
  229 
  230 ##############################
  231 # Test uniqind
  232 $x = pdl([0,1,2,2,0,1]);
  233 $y = $x->uniqind;
  234 eval { $c = all($y==pdl([0,1,3])) };
  235 is $@, '';
  236 ok $c, "uniqind" or diag "got: $y";
  237 is $y->ndims, 1, "uniqind";
  238 
  239 $y = pdl(1,1,1,1,1)->uniqind;         # SF bug 3076570
  240 ok(! $y->isempty);
  241 eval { $c = all($y==pdl([0])) };
  242 is $@, '';
  243 ok $c, "uniqind";
  244 is $y->ndims, 1, "uniqind, SF bug 3076570";
  245 
  246 ##############################
  247 # Test whereND
  248 $x = sequence(4,3,2);
  249 $y = pdl(0,1,1,0);
  250 $c = whereND($x,$y);
  251 ok(all(pdl($c->dims)==pdl(2,3,2))) and
  252 ok(all($c==pdl q[ [ [ 1  2] [ 5  6] [ 9 10] ]
  253                  [ [13 14] [17 18] [21 22] ] ]),
  254                                  "whereND [4]");
  255 $y = pdl q[ 0 0 1 1 ; 0 1 0 0 ; 1 0 0 0 ];
  256 $c = whereND($x,$y);
  257 ok(all(pdl($c->dims)==pdl(4,2))) and
  258 ok(all($c==pdl q[ 2  3  5  8 ; 14 15 17 20 ]),
  259                             "whereND [4,3]");
  260 $y = (random($x)<0.3);
  261 $c = whereND($x,$y);
  262 ok(all($c==where($x,$y)), "whereND vs where");
  263 # sf.net bug #3415115, whereND fails to handle all zero mask case
  264 $y = zeros(4);
  265 $c = whereND($x,$y);
  266 ok($c->isempty, 'whereND of all-zeros mask');
  267 # Make sure whereND functions as an lvalue:
  268 $x = sequence(4,3);
  269 $y = pdl(0, 1, 1, 1);
  270 eval { $x->whereND($y) *= -1 };
  271 is($@, '', 'using whereND in lvalue context does not croak');
  272 ok(all($x->slice("1:-1") < 0), 'whereND in lvalue context works');
  273 
  274 #Test fibonacci.
  275 my $fib=fibonacci(15);
  276 my $fib_ans = pdl(1,1,2,3,5,8,13,21,34,55,89,144,233,377,610);
  277 ok(all($fib == $fib_ans), 'Fibonacci sequence');
  278 
  279 #Test which_both.
  280 my $which_both_test=pdl(1,4,-2,0,5,0,1);
  281 my ($nonzero,$zero)=which_both($which_both_test);
  282 ok(all($nonzero==pdl(0,1,2,4,6)), 'Which_both nonzero indices');
  283 ok(all($zero==pdl(3,5)), 'Which_both zero indices');
  284 
  285 ###### Testing Begins #########
  286 my $im = new PDL [
  287   [ 1, 2,  3,  3 , 5],
  288   [ 2,  3,  4,  5,  6],
  289   [13, 13, 13, 13, 13],
  290   [ 1,  3,  1,  3,  1],
  291   [10, 10,  2,  2,  2,]
  292  ];
  293 my @minMax = $im->minmax;
  294 ok($minMax[0] == 1, "minmax min" );
  295 ok($minMax[1] == 13, "minmax max" );
  296 ok(($im x $im)->sum == 3429, "matrix multiplication" );
  297 my @statsRes = $im->stats;
  298 ok(tapprox($statsRes[0],5.36), "stats: mean" );
  299 ok(tapprox($statsRes[1],4.554), "stats: prms");
  300 ok(tapprox($statsRes[2],3), "stats: median");
  301 ok(tapprox($statsRes[3],1), "stats: min");
  302 ok(tapprox($statsRes[4],13), "stats: max");
  303 ok(tapprox($statsRes[6],4.462), "stats: rms");
  304 
  305 @statsRes = $im->short->stats; # Make sure that stats are promoted to floating-point
  306 ok(tapprox($statsRes[0],5.36), "stats: float mean");
  307 ok(tapprox($statsRes[1],4.554), "stats: float prms");
  308 ok(tapprox($statsRes[2],3), "stats: float median");
  309 ok(tapprox($statsRes[3],1), "stats: float min");
  310 ok(tapprox($statsRes[4],13), "stats: float max");
  311 ok(tapprox($statsRes[6],4.462), "stats: float rms");
  312 
  313 my $ones = ones(5,5);
  314 @statsRes = $im->stats($ones);
  315 ok(tapprox($statsRes[0],5.36), "stats: trivial weights mean" );
  316 ok(tapprox($statsRes[1],4.554), "stats: trivial weights prms" );
  317 ok(tapprox($statsRes[2],3), "stats: trivial weights median" );
  318 ok(tapprox($statsRes[3],1), "stats: trivial weights min" );
  319 ok(tapprox($statsRes[4],13), "stats: trivial weights max" );
  320 ok(tapprox($statsRes[6],4.462), "stats: trivial weights rms");
  321 
  322 # complex matmult
  323 my $cm1 = pdl('1 1+i 1');
  324 my $cm2 = pdl('2 3 i')->transpose;
  325 my $got = $cm1 x $cm2;
  326 ok all(approx $got, pdl('5+4i')), 'complex matmult' or diag $got;
  327 throws_ok { scalar $cm1->transpose x $cm2 } qr/mismatch/, 'good error on mismatch matmult';
  328 
  329 {
  330 my $pa = pdl [[ 1,  2,  3,  0],
  331       [ 1, -1,  2,  7],
  332       [ 1,  0,  0,  1]];
  333 my $pb = pdl [[1, 1],
  334      [0, 2],
  335      [0, 2],
  336      [1, 1]];
  337 my $pc = pdl [[ 1, 11],
  338       [ 8, 10],
  339       [ 2,  2]];
  340 my $res = $pa x $pb;
  341 ok(all approx($pc,$res)) or diag "got: $res";
  342 $res = null;
  343 matmult($pa, $pb, $res);
  344 ok(all(approx $pc,$res), 'res=null') or diag "got: $res";
  345 my $pa_sliced = $pa->dummy(0, 3)->dummy(-1, 3)->make_physical->slice('(1),,,(1)');
  346 $res = $pa_sliced x $pb;
  347 ok(all approx($pc,$res)) or diag "got: $res";
  348 $res = zeroes(2, 3);
  349 matmult($pa, $pb, $res);
  350 ok(all(approx $pc,$res), 'res=zeroes') or diag "got: $res";
  351 $res = ones(2, 3);
  352 matmult($pa, $pb, $res);
  353 ok(all(approx $pc,$res), 'res=ones') or diag "got: $res";
  354 my $eq = float [[1,1,1,1]];  # a 4,1-matrix ( 1 1 1 1 )
  355 # Check collapse: output should be a 1x2...
  356 ok(all approx($eq x $pb  , pdl([[2,6]]) )); # ([4x1] x [2x4] -> [1x2])
  357 # Check dimensional exception: mismatched dims should throw an error
  358 dies_ok {
  359 	my $pz = $pb x $eq; # [2x4] x [4x1] --> error (2 != 1)
  360 };
  361 {
  362 # Check automatic scalar multiplication
  363 my $pz;
  364 lives_ok { $pz = $pb x 2; };
  365 ok( all approx($pz,$pb * 2));
  366 }
  367 {
  368 my $pz;
  369 lives_ok { $pz = pdl(3) x $pb; };
  370 ok( all approx($pz,$pb * 3));
  371 }
  372 }
  373 
  374 # which ND test
  375 my $a1 = PDL->sequence(10,10,3,4);  
  376 
  377 ($x, $y, $z, my $w) = whichND($a1 == 203)->mv(0,-1)->dog;
  378 ok($a1->at($x->list,$y->list,$z->list,$w->list) == 203, "whichND" );
  379 
  380 $a1 = pdl(1,2,3,4);
  381 my $b1 = append($a1,2);
  382 ok(int(sum($b1))==12, "append");
  383 $b1 = append(null, null);
  384 ok !$b1->isnull, 'append(null, null) returns non-null';
  385 ok $b1->isempty, 'append(null, null) returns an empty';
  386 append(null, null, $b1);
  387 ok !$b1->isnull, 'append(null, null, b1) sets non-null';
  388 ok $b1->isempty, 'append(null, null, b1) sets an empty';
  389 $b1 = indx(1,2);
  390 is $b1->type, 'indx', '$indx_pdl is an indx pdl';
  391 $b1 = $b1->append(-1);
  392 is $b1->type, 'indx', 'append($indx_pdl, -1) returns an indx pdl';
  393 is $b1.'', '[1 2 -1]', 'append($indx_pdl, -1) correct content';
  394 
  395 # clip tests
  396 ok(tapprox($im->hclip(5)->sum,83), "hclip" );
  397 ok(tapprox($im->lclip(5)->sum,176), "lclip" );
  398 ok(tapprox($im->clip(5,7)->sum,140), "clip" );
  399 # with NaN badvalue
  400 $im = sequence(3);
  401 $im->badvalue(nan());
  402 $im->badflag(1);
  403 $im->set(1, nan());
  404 my $clipped = $im->lclip(0);
  405 is $clipped.'', '[0 BAD 2]', 'ISBAD() works when badvalue is NaN';
  406 
  407 # indadd Test:
  408 $a1 = pdl( 1,2,3);
  409 my $ind = pdl( 1,4,6);
  410 my $sum = zeroes(10);
  411 indadd($a1,$ind, $sum);
  412 ok(tapprox($sum->sum,6), "indadd" );
  413 
  414 #one2nd test
  415 $a1 = zeroes(3,4,5);
  416 my $indicies = pdl(0,1,4,6,23,58,59);
  417 ($x,$y,$z)=$a1->one2nd($indicies);
  418 ok(all( $x==pdl(0,1,1,0,2,1,2) ), "one2nd x");
  419 ok(all( $y==pdl(0,0,1,2,3,3,3) ), "one2nd y");
  420 ok(all( $z==pdl(0,0,0,0,1,4,4) ), "one2nd z");
  421 
  422 {
  423 my $yvalues =  (new PDL( 0..5))   - 20;
  424 my $xvalues = -(new PDL (0..5))*.5;
  425 my $x = new PDL(-2);
  426 is( $x->interpol($xvalues,$yvalues), -16 );
  427 }
  428 
  429 # Some of these tests are based upon those in Chapter 5 of Programming
  430 # Pearls, by J. Bentley
  431 {
  432 # choose a non-factor of two odd number for the length
  433 my $N = 723;
  434 
  435 my $ones = ones( $N );
  436 my $idx  = sequence( $N );
  437 my $x    = $idx * 10;
  438 
  439 # create ordered duplicates so can test insertion points. This creates
  440 # 7 sequential duplicates of the values 0-99
  441 my $ndup = 7;
  442 my $xdup = double long sequence( $ndup * 100 ) / $ndup;
  443 
  444 # get insertion points and values
  445 my ( $xdup_idx_insert_left, $xdup_idx_insert_right, $xdup_values ) = do {
  446 
  447     my ( $counts, $values ) = do { my @q = $xdup->rle; where( @q, $q[0] > 0 ) };
  448 
  449     ( $counts->cumusumover - $counts->at( 0 ), $counts->cumusumover, $values );
  450 
  451 };
  452 
  453 # The tests are table driven, with appropriate inputs and outputs for
  454 # forward and reverse sorted arrays.  The tests sort the input array
  455 # against itself, so we have a very good idea of which indices should
  456 # be returned.  Most of the tests use that.  There are also specific
  457 # tests for the endpoints as specified in the documentation, which
  458 # may be easier for humans to parse and validate.
  459 
  460 my %search = (
  461 
  462     sample => {
  463 
  464         all_the_same_element => $N - 1,    # finds right-most element
  465 
  466         forward => {
  467             idx      => $idx,
  468             x        => $x,
  469             equal    => $idx,
  470             nequal_m => $idx,
  471             nequal_p =>
  472               do { my $t = $idx + 1; $t->set( -1, $t->at( -1 ) - 1 ); $t },
  473             xdup => {
  474                 set    => $xdup,
  475                 idx    => $xdup_idx_insert_left,
  476                 values => $xdup_values,
  477             },
  478             docs => [
  479                 '          V <= xs[0] : i = 0                      ' => [ (  0, -1, 0 ),
  480 									  (  0,  0, 0 ),
  481 									],
  482                 'xs[0]  < V <= xs[-1] : i s.t. xs[i-1] < V <= xs[i]' => [ (  0,  1,  1 ),
  483 									  (  1,  0,  1 ),
  484 									  ( -1,  0, $N-1 ),
  485 									],
  486                 'xs[-1] < V           : i = $xs->nelem -1          ' => [ ( -1,  0, $N-1 ),
  487 									  ( -1,  1, $N-1 ),
  488 									],
  489             ],
  490         },
  491 
  492         reverse => {
  493             idx      => $idx,
  494             x        => $x->mslice( [ -1, 0 ] ),
  495             equal    => $idx,
  496             nequal_m => $idx,
  497             nequal_p => do { my $t = $idx - 1; $t->set( 0, 0 ); $t },
  498             xdup     => {
  499                 set => $xdup->slice( [ -1, 0 ] ),
  500                 idx => $xdup->nelem - 1 - $xdup_idx_insert_left,
  501                 values => $xdup_values,
  502             },
  503             docs => [
  504 		'          V > xs[0]  : i = 0                      ' => [(0,  1, 0) ],
  505 		'xs[0]  >= V > xs[-1] : i s.t. xs[i] >= V > xs[i+1]' => [(0,  0, 0),
  506 									 (0, -1, 0),
  507 									 (1,  0, 1),
  508 									],
  509 		'xs[-1] >= V          : i = $xs->nelem - 1         ' => [(-1,  0, $N-1),
  510 									 (-1, -1, $N-1),
  511 									],
  512             ],
  513        }
  514 
  515     },
  516 
  517     insert_leftmost => {
  518 
  519         all_the_same_element => 0,
  520 
  521         forward => {
  522             idx      => $idx,
  523             x        => $x,
  524             equal    => $idx,
  525             nequal_m => $idx,
  526             nequal_p => $idx + 1,
  527             xdup     => {
  528                 set    => $xdup,
  529                 idx    => $xdup_idx_insert_left,
  530                 values => $xdup_values,
  531             },
  532 	    docs => [
  533 		'         V <= xs[0]  : i = 0                      ' => [ ( 0, -1, 0 ),
  534 									  ( 0,  0,  0)
  535 									],
  536 		'xs[0]  < V <= xs[-1] : i s.t. xs[i-1] < V <= xs[i]' => [ ( 0, 1, 1 ),
  537 									  ( 1, 0, 1 ),
  538 									  ( -1, 0, $N-1 ),
  539 									],
  540 		'xs[-1] < V           : i = $xs->nelem             ' => [
  541 									 ( -1, 1, $N ),
  542 									],
  543 
  544 	    ],
  545         },
  546 
  547         reverse => {
  548             idx      => $idx,
  549             x        => $x->mslice( [ -1, 0 ] ),
  550             equal    => $idx,
  551             nequal_m => $idx,
  552             nequal_p => $idx - 1,
  553             xdup     => {
  554                 set => $xdup->mslice( [ -1, 0 ] ),
  555                 idx => $xdup->nelem - 1 - $xdup_idx_insert_left,
  556                 values => $xdup_values,
  557             },
  558            docs => [
  559 	       '          V >  xs[0]  : i = -1                     ' => [ ( 0,   1, -1 ), ],
  560 	       'xs[0]  >= V >= xs[-1] : i s.t. xs[i] >= V > xs[i+1]' => [ ( 0,   0,  0 ),
  561 									  ( 0,  -1,  0 ),
  562 									],
  563 	       'xs[-1] >= V           : i = $xs->nelem -1          ' => [ ( -1,  0, $N-1 ),
  564 									  ( -1, -1, $N-1 ),
  565 									],
  566 
  567            ],
  568         },
  569     },
  570 
  571     insert_rightmost => {
  572 
  573         all_the_same_element => $N,
  574 
  575         forward => {
  576             idx      => $idx,
  577             x        => $x,
  578             equal    => $idx + 1,
  579             nequal_m => $idx,
  580             nequal_p => $idx + 1,
  581             xdup     => {
  582                 set    => $xdup,
  583                 idx    => $xdup_idx_insert_right,
  584                 values => $xdup_values,
  585                 idx_offset => -1,   # returns index of element *after* the value
  586             },
  587 	    docs => [
  588 		'          V < xs[0]  : i = 0                      ' => [ ( 0, -1, 0 ) ],
  589 		'xs[0]  <= V < xs[-1] : i s.t. xs[i-1] <= V < xs[i]' => [ ( 0, 0, 1 ),
  590 									  ( 0, 1, 1 ),
  591 									  ( 1, 0, 2 ),
  592 									],
  593 		'xs[-1] <= V          : i = $xs->nelem             ' => [ ( -1, 0, $N ),
  594 									  ( -1, 1, $N ),
  595 									],
  596             ],
  597         },
  598 
  599         reverse => {
  600             idx      => $idx,
  601             x        => $x->mslice( [ -1, 0 ] ),
  602             equal    => $idx - 1,
  603             nequal_m => $idx,
  604             nequal_p => $idx - 1,
  605             xdup     => {
  606                 set => $xdup->mslice( [ -1, 0 ] ),
  607                 idx => $xdup->nelem - 1 - $xdup_idx_insert_right,
  608                 values => $xdup_values,
  609                 idx_offset => +1,   # returns index of element *after* the value
  610             },
  611 	    docs => [
  612 		'         V >= xs[0]  : i = -1                     ' => [ ( 0,   1, -1 ),
  613 									  ( 0,   0, -1 ),
  614 									],
  615 		'xs[0]  > V >= xs[-1] : i s.t. xs[i] >= V > xs[i+1]' => [ ( 0,  -1, 0 ),
  616 									  ( -1,  1, $N-2 ),
  617 									  ( -1,  0, $N-2 ),
  618 									],
  619 		'xs[-1] > V           : i = $xs->nelem -1          ' => [ ( -1,  -1, $N-1 ) ]
  620             ],
  621         },
  622     },
  623 
  624     match => {
  625 
  626         all_the_same_element => ( $N ) >> 1,
  627 
  628         forward => {
  629             idx      => $idx,
  630             x        => $x,
  631             equal    => $idx,
  632             nequal_m => -( $idx + 1 ),
  633             nequal_p => -( $idx + 1 + 1 ),
  634             xdup     => {
  635                 set    => $xdup,
  636                 values => $xdup_values,
  637             },
  638 	    docs => [
  639 		'V < xs[0]  : i = -1' => [ ( 0,   -1, -1 ), ],
  640 		'V == xs[n] : i = n' => [ ( 0,  0, 0 ),
  641 					  ( -1, 0, $N-1 ) ],
  642 		'xs[0] > V > xs[-1], V != xs[n] : -(i+1) s.t. xs[i] > V > xs[i+1]' => [ ( 0,   1, -( 1 + 1)  ),
  643 											( 1,  -1, -( 1 + 1 ) ),
  644 											( 1,   1, -( 2 + 1 ) ),
  645 											( -1,  -1, -( $N - 1 + 1 ) ),
  646 										      ],
  647 		' V > xs[-1] : -($xs->nelem - 1 + 1)' => [ ( -1,   1, -( $N + 1)  ), ]
  648             ],
  649         },
  650 
  651         reverse => {
  652             idx      => $idx,
  653             x        => $x->mslice( [ -1, 0 ] ),
  654             equal    => $idx,
  655             nequal_m => -( $idx + 1 ),
  656             nequal_p => -( $idx + 1 - 1 ),
  657             xdup     => {
  658                 set => $xdup->mslice( [ -1, 0 ] ),
  659                 values => $xdup_values,
  660             },
  661 	    docs => [
  662 		'V > xs[0]  : i = 0' => [ ( 0,  1, 0 ), ],
  663 		'V == xs[n] : i = n' => [ ( 0,  0, 0 ),
  664 					  ( -1, 0, $N-1 ) ],
  665 		'xs[0] < V < xs[-1], V != xs[n] : -(i+1) s.t. xs[i-1] > V > xs[i]' => [ ( 0,  -1, -( 0 + 1)  ),
  666 											( 1,   1, -( 0 + 1 ) ),
  667 											( 1,  -1, -( 1 + 1 ) ),
  668 											( -1,  -1, -( $N - 1 + 1 ) ),
  669 										      ],
  670 		' xs[-1] > V: -($xs->nelem - 1 + 1)' => [ ( -1,   -1, -( $N - 1 + 1)  ), ]
  671             ],
  672         },
  673     },
  674 
  675     bin_inclusive => {
  676 
  677         all_the_same_element => $N - 1,
  678 
  679         forward => {
  680             idx      => $idx,
  681             x        => $x,
  682             equal    => $idx,
  683             nequal_m => $idx - 1,
  684             nequal_p => $idx,
  685             xdup     => {
  686                 set    => $xdup,
  687                 idx    => $xdup_idx_insert_left + $ndup - 1,
  688                 values => $xdup_values,
  689             },
  690 	    docs => [
  691 		'          V < xs[0]  : i = -1                     ' => [ ( 0, -1, -1 ), ],
  692 		'xs[0]  <= V < xs[-1] : i s.t. xs[i] <= V < xs[i+1]' => [ ( 0,  0,  0 ),
  693 									  ( 0,  1,  0 ),
  694 									  ( 1, -1,  0 ),
  695 									  ( 1,  0,  1 ),
  696 									  ( -1, -1, $N-2 ),
  697 									],
  698 		'xs[-1] <= V          : i = $xs->nelem - 1         ' => [
  699 									  ( -1, 0,  $N-1 ),
  700 									  ( -1, 1,  $N-1 ),
  701 									]
  702             ],
  703         },
  704 
  705         reverse => {
  706             idx      => $idx,
  707             x        => $x->mslice( [ -1, 0 ] ),
  708             equal    => $idx,
  709             nequal_m => $idx + 1,
  710             nequal_p => $idx,
  711             xdup     => {
  712                 set => $xdup->mslice( [ -1, 0 ] ),
  713                 idx => $xdup->nelem - ( 1 + $xdup_idx_insert_left + $ndup - 1 ),
  714                 values => $xdup_values,
  715             },
  716 	    docs => [
  717 		'          V >= xs[0]  : i = 0                        ' => [ (0, 1, 0 ),
  718 									     (0, 0, 0 )
  719 									 ],
  720 		'xs[0]  >  V >= xs[-1] : i s.t. xs[i+1] > V >= xs[i]' => [ ( 0, -1, 1 ),
  721 									   ( 1,  1, 1 ),
  722 									   ( 1,  0, 1 ),
  723 									   ( 1, -1, 2 ),
  724 									   ( -1, 0, $N-1 ),
  725 									 ],
  726 		'xs[-1] >  V           : i = $xs->nelem -1          ' => [ ( -1, -1, $N ) ],
  727             ],
  728         },
  729     },
  730 
  731     bin_exclusive => {
  732 
  733         all_the_same_element => -1,
  734 
  735         forward => {
  736             idx      => $idx,
  737             x        => $x,
  738             equal    => $idx - 1,
  739             nequal_m => $idx - 1,
  740             nequal_p => $idx,
  741             xdup     => {
  742                 set        => $xdup,
  743                 idx        => $xdup_idx_insert_left - 1,
  744                 values     => $xdup_values,
  745                 idx_offset => 1,
  746             },
  747 	    docs => [
  748 		'          V <= xs[0]  : i = -1                     ' => [ ( 0, -1, -1 ),
  749 									   ( 0,  0, -1 ),
  750 									 ],
  751 		'xs[0]  <  V <= xs[-1] : i s.t. xs[i] < V <= xs[i+1]' => [ ( 0,  1, 0 ),
  752 									   ( 1, -1, 0 ),
  753 									   ( 1,  0, 0 ),
  754 									   ( 1,  1, 1 ),
  755 									   ( -1, -1, $N-2 ),
  756 									   ( -1, 0, $N-2 ),
  757 									],
  758 		'xs[-1] <  V           : i = $xs->nelem - 1         ' => [
  759 									  ( -1, 1, $N-1 ),
  760 									 ],
  761             ],
  762         },
  763 
  764         reverse => {
  765             idx      => $idx,
  766             x        => $x->mslice( [ -1, 0 ] ),
  767             equal    => $idx + 1,
  768             nequal_m => $idx + 1,
  769             nequal_p => $idx,
  770             xdup     => {
  771                 set => $xdup->mslice( [ -1, 0 ] ),
  772                 idx => $xdup->nelem - ( 1 + $xdup_idx_insert_left - 1 ),
  773                 values     => $xdup_values,
  774                 idx_offset => -1,
  775             },
  776 	    docs => [
  777 		'          V >  xs[0]  : i = 0                      ' => [ ( 0,  1, 0 ), ],
  778 		'xs[0]  >  V >  xs[-1] : i s.t. xs[i-1] >= V > xs[i]' => [ ( 0,  0, 1 ),
  779 									   ( 0, -1, 1 ),
  780 									   ( -1, 1, $N-1 ),
  781 									 ],
  782 		'xs[-1] >= V           : i = $xs->nelem -1          ' => [ ( -1, 0, $N ),
  783 									   ( -1, -1, $N ),
  784 									 ],
  785 	    ],
  786         },
  787     },
  788 
  789 );
  790 
  791 for my $mode (
  792     sort keys %search
  793   )
  794 {
  795 
  796     my $data   = $search{$mode};
  797 
  798     subtest $mode => sub {
  799 
  800         my ( $got, $exp );
  801         for my $sort_direction ( qw[ forward reverse ] ) {
  802 
  803             subtest $sort_direction => sub {
  804 
  805 		my $so = $data->{$sort_direction}
  806 		  or plan( skip_all => "not testing $sort_direction!\n" );
  807 
  808                 ok(
  809                     all(
  810                         ( $got = vsearch( $so->{x}, $so->{x}, { mode => $mode } ) )
  811 			==
  812 			( $exp = $so->{equal} )
  813                     ),
  814                     'equal elements'
  815                 ) or diag "got     : $got\nexpected: $exp\n";
  816 
  817                 ok(
  818                     all(
  819                         ( $got = vsearch( $so->{x} - 5, $so->{x}, { mode => $mode } ) )
  820                         ==
  821 			( $exp = $so->{nequal_m} )
  822                     ),
  823                     'non-equal elements x[i] < xs[i] (check lower bound)'
  824                 ) or diag "got     : $got\nexpected: $exp\n";
  825 
  826                 ok(
  827                     all(
  828                         ( $got = vsearch( $so->{x} + 5, $so->{x}, { mode => $mode } ) )
  829                         ==
  830 			( $exp = $so->{nequal_p} )
  831                     ),
  832                     'non-equal elements x[i] > xs[i] (check upper bound)'
  833                 ) or diag "got     : $got\nexpected: $exp\n";
  834 
  835 
  836 		# duplicate testing.
  837 
  838 		# check for values. note that the rightmost routine returns
  839 		# the index of the element *after* the last duplicate
  840 		# value, so we need an offset
  841 		ok(
  842 		    all(
  843 			( $got = $so->{xdup}{set}->index( vsearch( $so->{xdup}{values}, $so->{xdup}{set}, { mode => $mode } )
  844 							                 + ($so->{xdup}{idx_offset} || 0) ) )
  845 			==
  846 			( $exp = $so->{xdup}{values} )
  847 		    ),
  848 		    'duplicates values'
  849 		) or diag "got     : $got\nexpected: $exp\n";
  850 
  851 		# if there are guarantees about which duplicates are returned, test it
  852 		if ( exists $so->{xdup}{idx} ) {
  853 
  854 		    ok(
  855 			all(
  856 			    ( $got = vsearch( $so->{xdup}{values}, $so->{xdup}{set}, { mode => $mode } ) )
  857 			    ==
  858 			    ( $exp = $so->{xdup}{idx} )
  859 			),
  860 			'duplicate indices'
  861 		    ) or diag "got     : $got\nexpected: $exp\n";
  862 
  863 		}
  864 
  865 		if ( exists $so->{docs} ) {
  866 
  867 		    while( my ($label, $inputs ) = splice( @{$so->{docs}}, 0, 2 )  ) {
  868 
  869 			while( @$inputs ) {
  870 
  871 			    my ( $idx, $offset, $exp ) = splice( @$inputs, 0, 3 );
  872 			    my $value = $so->{x}->at($idx) + $offset;
  873 
  874 			    is ( $got = ( vsearch( $value, $so->{x}, { mode => $mode } )->sclr), $exp, "$label: ($idx, $offset)" );
  875 
  876 			}
  877 		    }
  878 		}
  879 
  880 
  881             };
  882         }
  883 
  884         ok(
  885             all(
  886                 ( $got = vsearch( $ones, $ones, { mode => $mode } ) )
  887                 ==
  888                 ( $exp = $data->{all_the_same_element} )
  889             ),
  890             'all the same element'
  891         ) or diag "got     : $got\nexpected: $exp\n";
  892     };
  893 
  894 }
  895 
  896 # test vsearch API to ensure backwards compatibility
  897 {
  898     my $vals = random( 100 );
  899     my $xs = sequence(100) / 99;
  900 
  901     # implicit output ndarray
  902     my $indx0 = vsearch( $vals, $xs );
  903 
  904     my $ret = vsearch( $vals, $xs, my $indx1 = PDL->null() );
  905 
  906     is( $ret, undef, "no return from explicit output ndarray" );
  907 
  908     ok ( all ( $indx0 == $indx1 ),
  909 	 'explicit ndarray == implicit ndarray' );
  910 }
  911 }
  912 
  913 my $vdim = 4;
  914 my $v1 = zeroes($vdim);
  915 my $v2 = pdl($v1);
  916 $v2->set(-1,1);
  917 
  918 ok $v1->cmpvec($v2)<0, "cmpvec:1d:<";
  919 ok $v2->cmpvec($v1)>0, "cmpvec:1d:>";
  920 is $v1->cmpvec($v1)->sclr, 0, "cmpvec:1d:==";
  921 
  922 ##-- 4..5: qsortvec, qsortveci
  923 my $p2d  = pdl([[1,2],[3,4],[1,3],[1,2],[3,3]]);
  924 
  925 ok all(approx($p2d->qsortvec, pdl(long,[[1,2],[1,2],[1,3],[3,3],[3,4]]))), "qsortvec";
  926 ok all(approx($p2d->dice_axis(1,$p2d->qsortveci), $p2d->qsortvec)), "qsortveci";
  927 
  928 my $which = pdl(long,[[0,0],[0,0],[0,1],[0,1],[1,0],[1,0],[1,1],[1,1]]);
  929 my $find  = $which->slice(",0:-1:2");
  930 
  931 ok all(approx($find->vsearchvec($which), pdl(long,[0,2,4,6]))), "vsearchvec():match";
  932 ok all(pdl([-1,-1])->vsearchvec($which)==0), "vsearchvec():<<";
  933 ok all(pdl([2,2])->vsearchvec($which)==$which->dim(1)-1), "vsearchvec():>>";
  934 
  935 my $vtype = long;
  936 my $universe = pdl($vtype,[ [0,0],[0,1],[1,0],[1,1] ]);
  937 $v1 = $universe->dice_axis(1,pdl([0,1,2]));
  938 $v2 = $universe->dice_axis(1,pdl([1,2,3]));
  939 
  940 ($c,my $nc) = $v1->unionvec($v2);
  941 ok all(approx($c, pdl($vtype, [ [0,0],[0,1],[1,0],[1,1],[0,0],[0,0] ]))), "unionvec:list:c";
  942 is $nc, $universe->dim(1), "unionvec:list:nc";
  943 my $cc = $v1->unionvec($v2);
  944 ok all(approx($cc, $universe)), "unionvec:scalar";
  945 
  946 ($c,$nc) = $v1->intersectvec($v2);
  947 ok all(approx($c, pdl($vtype, [ [0,1],[1,0],[0,0] ]))), "intersectvec:list:c";
  948 is $nc->sclr, 2, "intersectvec:list:nc";
  949 $cc = $v1->intersectvec($v2);
  950 ok all(approx($cc, $universe->slice(",1:2"))), "intersectvec:scalar";
  951 
  952 ($c,$nc) = $v1->setdiffvec($v2);
  953 ok all(approx($c, pdl($vtype, [ [0,0], [0,0],[0,0] ]))), "setdiffvec:list:c";
  954 is $nc, 1, "setdiffvec:list:nc";
  955 $cc = $v1->setdiffvec($v2);
  956 ok all(approx($cc, pdl($vtype, [[0,0]]))), "setdiffvec:scalar";
  957 
  958 my $all = sequence(20);
  959 my $amask = ($all % 2)==0;
  960 my $bmask = ($all % 3)==0;
  961 my $a   = $all->where($amask);
  962 my $b   = $all->where($bmask);
  963 
  964 ok all(approx(scalar($a->union_sorted($b)), $all->where($amask | $bmask))), "union_sorted";
  965 ok all(approx(scalar($a->intersect_sorted($b)),  $all->where($amask & $bmask))), "intersect_sorted";
  966 ok all(approx(scalar($a->setdiff_sorted($b)), $all->where($amask & $bmask->not))), "setdiff_sorted";
  967 
  968 ##--------------------------------------------------------------
  969 ## dim-checks and implicit broadcast dimensions
  970 ##  + see https://github.com/moocow-the-bovine/PDL-VectorValued/issues/4
  971 
  972 sub test_broadcast_dimensions {
  973   ##-- unionvec
  974   my $empty = zeroes(3,0);
  975   my $uw = pdl([[-3,-2,-1],[1,2,3]]);
  976   my $wx = pdl([[1,2,3],[4,5,6]]);
  977   my $xy = pdl([[4,5,6],[7,8,9]]);
  978 
  979   # unionvec: basic
  980   ok all(approx(scalar($uw->unionvec($wx)), pdl([[-3,-2,-1],[1,2,3],[4,5,6]]))), "unionvec - broadcast dims - uw+wx";
  981   ok all(approx(scalar($uw->unionvec($xy)), pdl([[-3,-2,-1],[1,2,3],[4,5,6],[7,8,9]]))), "unionvec - broadcast dims - uw+xy";
  982   ok all(approx(scalar($empty->unionvec($wx)), $wx)), "unionvec - broadcast dims - 0+wx";
  983   ok all(approx(scalar($wx->unionvec($empty)), $wx)), "unionvec - broadcast dims - wx+0";
  984   ok all(approx(scalar($empty->unionvec($empty)), $empty)), "unionvec - broadcast dims - 0+0";
  985 
  986   # unionvec: broadcasting
  987   my $k = 2;
  988   my $kempty = $empty->slice(",,*$k");
  989   my $kuw = $uw->slice(",,*$k");
  990   my $kwx = $wx->slice(",,*$k");
  991   my $kxy = $xy->slice(",,*$k");
  992   ok all(approx(scalar($kuw->unionvec($wx)), pdl([[-3,-2,-1],[1,2,3],[4,5,6]])->slice(",,*$k"))), "unionvec - broadcast dims - uw(*k)+wx";
  993   ok all(approx(scalar($kuw->unionvec($xy)), pdl([[-3,-2,-1],[1,2,3],[4,5,6],[7,8,9]])->slice(",,*$k"))), "unionvec - broadcast dims - uw(*k)+xy";
  994   ok all(approx(scalar($kempty->unionvec($wx)), $kwx)), "unionvec - broadcast dims - 0(*k)+wx";
  995   ok all(approx(scalar($kwx->unionvec($empty)), $kwx)), "unionvec - broadcast dims - wx(*k)+0";
  996   ok all(approx(scalar($kempty->unionvec($empty)), $kempty)), "unionvec - broadcast dims - 0(*k)+0";
  997 
  998   ##-- intersectvec
  999 
 1000   my $needle0 = pdl([[-3,-2,-1]]);
 1001   my $needle1 = pdl([[1,2,3]]);
 1002   my $needles = pdl([[-3,-2,-1],[1,2,3]]);
 1003   my $haystack = pdl([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]);
 1004 
 1005   # intersectvec: basic
 1006   ok all(approx(scalar($needle0->intersectvec($haystack)), $empty)), "intersectvec - broadcast dims - needle0&haystack";
 1007   ok all(approx(scalar($needle1->intersectvec($haystack)), $needle1)), "intersectvec - broadcast dims - needle1&haystack";
 1008   ok all(approx(scalar($needles->intersectvec($haystack)), $needle1)), "intersectvec - broadcast dims - needles&haystack";
 1009   ok all(approx(scalar($haystack->intersectvec($haystack)), $haystack)), "intersectvec - broadcast dims - haystack&haystack";
 1010   ok all(approx(scalar($haystack->intersectvec($empty)), $empty)), "intersectvec - broadcast dims - haystack&empty";
 1011   ok all(approx(scalar($empty->intersectvec($haystack)), $empty)), "intersectvec - broadcast dims - empty&haystack";
 1012 
 1013   # intersectvec: broadcasting
 1014   my $kneedle0 = $needle0->slice(",,*$k");
 1015   my $kneedle1 = $needle1->slice(",,*$k");
 1016   my $kneedles = pdl([[[-3,-2,-1]],[[1,2,3]]]);
 1017   my $khaystack = $haystack->slice(",,*$k");
 1018   ok all(approx(scalar($kneedle0->intersectvec($haystack)), $kempty)), "intersectvec - broadcast dims - needle0(*k)&haystack";
 1019   ok all(approx(scalar($kneedle1->intersectvec($haystack)), $kneedle1)), "intersectvec - broadcast dims - needle1(*k)&haystack";
 1020   ok all(approx(
 1021 	scalar($kneedles->intersectvec($haystack)),
 1022 	pdl([[[0,0,0]],[[1,2,3]]]))), "intersectvec - broadcast dims - needles(*k)&haystack";
 1023   ok all(approx(scalar($khaystack->intersectvec($haystack)), $khaystack)), "intersectvec - broadcast dims - haystack(*k)&haystack";
 1024   ok all(approx(scalar($khaystack->intersectvec($empty)), $kempty)), "intersectvec - broadcast dims - haystack(*k)&empty";
 1025   ok all(approx(scalar($kempty->intersectvec($haystack)), $kempty)), "intersectvec - broadcast dims - empty(*k)&haystack";
 1026 
 1027   ##-- setdiffvec
 1028 
 1029   # setdiffvec: basic
 1030   ok all(approx(scalar($haystack->setdiffvec($needle0)), $haystack)), "setdiffvec - broadcast dims - haystack-needle0";
 1031   ok all(approx(scalar($haystack->setdiffvec($needle1)), $haystack->slice(",1:-1"))), "setdiffvec - broadcast dims - haystack-needle1";
 1032   ok all(approx(scalar($haystack->setdiffvec($needles)), $haystack->slice(",1:-1"))), "setdiffvec - broadcast dims - haystack-needles";
 1033   ok all(approx(scalar($haystack->setdiffvec($haystack)), $empty)), "setdiffvec - broadcast dims - haystack-haystack";
 1034   ok all(approx(scalar($haystack->setdiffvec($empty)), $haystack)), "setdiffvec - broadcast dims - haystack-empty";
 1035   ok all(approx(scalar($empty->setdiffvec($haystack)), $empty)), "setdiffvec - broadcast dims - empty-haystack";
 1036 
 1037   # setdiffvec: broadcasting
 1038   ok all(approx(scalar($khaystack->setdiffvec($needle0)), $khaystack)), "setdiffvec - broadcast dims - haystack(*k)-needle0";
 1039   ok all(approx(scalar($khaystack->setdiffvec($needle1)), $khaystack->slice(",1:-1,"))), "setdiffvec - broadcast dims - haystack(*k)-needle1";
 1040   ok all(approx(scalar($khaystack->setdiffvec($needles)), $khaystack->slice(",1:-1,"))), "setdiffvec - broadcast dims - haystack(*k)-needles";
 1041   ok all(approx(scalar($khaystack->setdiffvec($haystack)), $kempty)), "setdiffvec - broadcast dims - haystack(*k)-haystack";
 1042   ok all(approx(scalar($khaystack->setdiffvec($empty)), $khaystack)), "setdiffvec - broadcast dims - haystack(*k)-empty";
 1043   ok all(approx(scalar($kempty->setdiffvec($haystack)), $kempty)), "setdiffvec - broadcast dims - empty(*k)-haystack";
 1044 }
 1045 test_broadcast_dimensions();
 1046 
 1047 ## intersectvec tests as suggested by ETJ/mowhawk2
 1048 ##  + see https://github.com/moocow-the-bovine/PDL-VectorValued/issues/4
 1049 sub test_intersect_implicit_dims {
 1050   # intersectvec: from ETJ/mowhawk2 a la https://stackoverflow.com/a/71446817/3857002
 1051   my $toto = pdl([1,2,3], [4,5,6]);
 1052   my $titi = pdl(1,2,3);
 1053   my $notin = pdl(7,8,9);
 1054   my ($c);
 1055 
 1056   ok all(approx($c=intersectvec($titi,$toto), [[1,2,3]])), 'intersectvec - implicit dims - titi&toto';
 1057   ok all(approx($c=intersectvec($notin,$toto), zeroes(3,0))), 'intersectvec - implicit dims - notin&toto';
 1058   ok all(approx($c=intersectvec($titi->dummy(1), $toto), [[1,2,3]])), 'intersectvec - implicit dims - titi(*1)&toto';
 1059   ok all(approx($c=intersectvec($notin->dummy(1), $toto), zeroes(3,0))), 'intersectvec - implicit dims - notin(*1)&toto';
 1060 
 1061   my $needle0_in = pdl([1,2,3]); # 3
 1062   my $needle0_notin = pdl([9,9,9]); # 3
 1063   my $needle_in = $needle0_in->dummy(1);  # 3x1: [[1 2 3]]
 1064   my $needle_notin = $needle0_notin->dummy(1); # 3x1: [[-3 -2 -1]]
 1065   my $needles = pdl([[1,2,3],[9,9,9]]); # 3x2: $needle0_in->cat($needle0_notin)
 1066   my $haystack = pdl([[1,2,3],[4,5,6]]); # 3x2
 1067 
 1068   sub intersect_ok {
 1069     my ($label, $a,$b, $c_want,$nc_want,$c_sclr_want) = @_;
 1070     my ($c, $nc) = intersectvec($a,$b);
 1071     my $c_sclr = intersectvec($a,$b);
 1072     ok all(approx($c, $c_want)), "$label - result";
 1073     ok all(approx($nc, $nc_want)), "$label - counts";
 1074     ok all(approx($c_sclr, $c_sclr_want)), "$label - scalar";
 1075   }
 1076 
 1077   intersect_ok('intersectvec - implicit dims - needle0_in&haystack',
 1078 	       $needle0_in, $haystack,
 1079 	       [[1,2,3]], 1, [[1,2,3]]
 1080 	      );
 1081   intersect_ok('intersectvec - implicit dims - needle_in&haystack',
 1082 	       $needle_in, $haystack,
 1083 	       [[1,2,3]], 1, [[1,2,3]]
 1084 	      );
 1085 
 1086   intersect_ok('intersectvec - implicit dims - needle0_notin&haystack',
 1087 	       $needle0_notin, $haystack,
 1088 	       [[0,0,0]], 0, zeroes(3,0)
 1089 	      );
 1090   intersect_ok('intersectvec - implicit dims - needle_notin&haystack',
 1091 	       $needle_notin, $haystack,
 1092 	       [[0,0,0]], 0, zeroes(3,0)
 1093 	      );
 1094 
 1095   intersect_ok('intersectvec - implicit dims - needles&haystack',
 1096 	       $needles, $haystack,
 1097 	       [[1,2,3],[0,0,0]], 1, [[1,2,3]]
 1098 	      );
 1099 
 1100   # now we want to know whether each needle is "in" one by one, not really
 1101   # a normal intersect, so we insert a dummy in haystack in order to broadcast
 1102   # the "nc" needs to come back as a 4x2
 1103   my $needles8 = pdl( [[[1,2,3],[4,5,6],[8,8,8],[8,8,8]],
 1104 		       [[4,5,6],[9,9,9],[1,2,3],[9,9,9]]]); # 3x4x2
 1105 
 1106   # need to manipulate above into suitable inputs for intersect to get right output
 1107   # + dummy dim here also ensures singleton query-vector-sets are (trivially) sorted
 1108   my $needles8x = $needles8->slice(",*1,,"); # 3x*x4x2 # dummy of size 1 inserted in dim 1
 1109 
 1110   # haystack: no changes needed; don't need same number of dims, broadcast engine will add dummy/1s at top
 1111   my $haystack8 = $haystack;
 1112   my $c_want8 = [
 1113 		 [[[1,2,3]],[[4,5,6]],[[0,0,0]],[[0,0,0]]],
 1114 		 [[[4,5,6]],[[0,0,0]],[[1,2,3]],[[0,0,0]]],
 1115 		];
 1116   my $nc_want8 = [[1,1,0,0],
 1117 		  [1,0,1,0]];
 1118 
 1119   intersect_ok('intersectvec - implicit dims - needles8x&haystack8',
 1120 	       $needles8x, $haystack8,
 1121 	       $c_want8, $nc_want8, $c_want8
 1122 	      );
 1123 }
 1124 test_intersect_implicit_dims();
 1125 
 1126 ## dim-checks and implicit broadcast dimensions
 1127 ##  + analogous to https://github.com/moocow-the-bovine/PDL-VectorValued/issues/4
 1128 sub test_v_broadcast_dimensions {
 1129   # data: basic
 1130   my $empty = zeroes(0);
 1131   my $v1_2 = pdl([1,2]);
 1132   my $v3_4 = pdl([3,4]);
 1133   my $v1_4 = $v1_2->cat($v3_4)->flat;
 1134 
 1135   # data: broadcasting
 1136   my $k = 2;
 1137   my $kempty = $empty->slice(",*$k");
 1138   my $kv1_2 = $v1_2->slice(",*$k");
 1139   my $kv3_4 = $v3_4->slice(",*$k");
 1140   my $kv1_4 = $v1_4->slice(",*$k");
 1141 
 1142   #-- union_sorted
 1143   ok all(approx(scalar($v1_2->union_sorted($v3_4)), $v1_4)), "union_sorted - broadcast dims - 12+34";
 1144   ok all(approx(scalar($v3_4->union_sorted($v1_4)), $v1_4)), "union_sorted - broadcast dims - 34+1234";
 1145   ok all(approx(scalar($empty->union_sorted($v1_4)), $v1_4)), "union_sorted - broadcast dims - 0+1234";
 1146   ok all(approx(scalar($v1_4->union_sorted($empty)), $v1_4)), "union_sorted - broadcast dims - 1234+0";
 1147   ok all(approx(scalar($empty->union_sorted($empty)), $empty)), "union_sorted - broadcast dims - 0+0";
 1148   #
 1149   ok all(approx(scalar($kv1_2->union_sorted($v3_4)), $kv1_4)), "union_sorted - broadcast dims - 12(*k)+34";
 1150   ok all(approx(scalar($kv3_4->union_sorted($v1_4)), $kv1_4)), "union_sorted - broadcast dims - 34(*k)+1234";
 1151   ok all(approx(scalar($kempty->union_sorted($v1_4)), $kv1_4)), "union_sorted - broadcast dims - 0(*k)+1234";
 1152   ok all(approx(scalar($kv1_4->union_sorted($empty)), $kv1_4)), "union_sorted - broadcast dims - 1234(*k)+0";
 1153   ok all(approx(scalar($kempty->union_sorted($empty)), $kempty)), "union_sorted - broadcast dims - 0(*k)+0";
 1154 
 1155   #-- intersect_sorted
 1156   ok all(approx(scalar($v1_2->intersect_sorted($v3_4)), $empty)), "intersect_sorted - broadcast dims - 12&34";
 1157   ok all(approx(scalar($v3_4->intersect_sorted($v1_4)), $v3_4)), "intersect_sorted - broadcast dims - 34&1234";
 1158   ok all(approx(scalar($empty->intersect_sorted($v1_4)), $empty)), "intersect_sorted - broadcast dims - 0&1234";
 1159   ok all(approx(scalar($v1_4->intersect_sorted($empty)), $empty)), "intersect_sorted - broadcast dims - 1234&0";
 1160   ok all(approx(scalar($empty->intersect_sorted($empty)), $empty)), "intersect_sorted - broadcast dims - 0&0";
 1161   #
 1162   ok all(approx(scalar($kv1_2->intersect_sorted($v3_4)), $kempty)), "intersect_sorted - broadcast dims - 12(*k)&34";
 1163   ok all(approx(scalar($kv3_4->intersect_sorted($v1_4)), $kv3_4)), "intersect_sorted - broadcast dims - 34(*k)&1234";
 1164   ok all(approx(scalar($kempty->intersect_sorted($v1_4)), $kempty)), "intersect_sorted - broadcast dims - 0(*k)&1234";
 1165   ok all(approx(scalar($kv1_4->intersect_sorted($empty)), $kempty)), "intersect_sorted - broadcast dims - 1234(*k)&0";
 1166   ok all(approx(scalar($kempty->intersect_sorted($empty)), $kempty)), "intersect_sorted - broadcast dims - 0(*k)&0";
 1167 
 1168   #-- setdiff_sorted
 1169   ok all(approx(scalar($v1_2->setdiff_sorted($v3_4)), $v1_2)), "setdiff_sorted - broadcast dims - 12-34";
 1170   ok all(approx(scalar($v3_4->setdiff_sorted($v1_4)), $empty)), "setdiff_sorted - broadcast dims - 34-1234";
 1171   ok all(approx(scalar($v1_4->setdiff_sorted($empty)), $v1_4)), "setdiff_sorted - broadcast dims - 1234-0";
 1172   ok all(approx(scalar($empty->setdiff_sorted($v1_4)), $empty)), "setdiff_sorted - broadcast dims - 0-1234";
 1173   ok all(approx(scalar($empty->setdiff_sorted($empty)), $empty)), "setdiff_sorted - broadcast dims - 0-0";
 1174   #
 1175   ok all(approx(scalar($kv1_2->setdiff_sorted($v3_4)), $kv1_2)), "setdiff_sorted - broadcast dims - 12(*k)-34";
 1176   ok all(approx(scalar($kv3_4->setdiff_sorted($v1_4)), $kempty)), "setdiff_sorted - broadcast dims - 34(*k)-1234";
 1177   ok all(approx(scalar($kv1_4->setdiff_sorted($empty)), $kv1_4)), "setdiff_sorted - broadcast dims - 1234(*k)-0";
 1178   ok all(approx(scalar($kempty->setdiff_sorted($v1_4)), $kempty)), "setdiff_sorted - broadcast dims - 0(*k)-1234";
 1179   ok all(approx(scalar($kempty->setdiff_sorted($empty)), $kempty)), "setdiff_sorted - broadcast dims - 0(*k)-0";
 1180 }
 1181 test_v_broadcast_dimensions();
 1182 
 1183 done_testing;