vr your code is superior than the long code I have posted!

If I may add: using oddpctover() might be preferred because it does not interpolate when there is no data at the exact percentile position.

Regarding the time difference when running with and without "use threads", I have discovered that avg() is the culprit. If you use x-> index( $idx )->at(0) rather than x-> index( $idx )->avg the performance is the same (which means idx() is also excluded as possible cause).

Regarding the time difference when running with and without "use threads", I have discovered that avg() is the culprit.

Hm-m, it's not what I'm observing here. Setting $m = 30_000 and replacing avg with at(0), these are results of 3 runs without thread and in a thread:

14.9562268257141
14.9582891464233
14.8853561878204

11.9686307907104
12.0527169704437
12.0850310325623
[download]

And then replacing all 3 lines of loop block with just simple

my $y = $x-> index( sequence $n )-> at( 0 );
[download]

7.74871516227722
7.80155396461487
7.71721601486206

4.92977499961853
4.87044596672058
4.87968802452087
[download]

So I'd say it's something strange with index going on. I'm not appealing to anyone for investigation :), observable speed difference may depend very much on hardware, let it be a murky PDL mystery.

And another entertaining (but dangerous, maybe) bit of PDL trivia: above I sketched a parallelization solution with threads and random. Except, if there were more than 1 worker thread, it won't work as expected. The random documentation says Perl's srand can be used to seed, and one may assume that random relies on Perl's RNG, including it (automatically) seeds in each thread that's started. Consider:

use strict;
use warnings;
use feature 'say';
use threads;
use PDL;

PDL::no_clone_skip_warning;

srand;
async sub{ say random( 1 )-> at( 0 )}
    for 1 .. 5;
$_-> join for threads-> list;

say rand;

__END__

0.851411183904023
0.851411183904023
0.851411183904023
0.851411183904023
0.851411183904023
0.851411183904023
[download]

That's some randomness. Try to say rand instead. So, one must to explicitly call srand at the start of a thread, if using threads and PDL's random.

Hi, vr

Tonight came across your post and modified your demonstration to run with 4 threads.

# https://www.perlmonks.org/?node_id=1214227

use strict;
use warnings;
use feature 'say';

use PDL;
use PDL::Parallel::threads qw(retrieve_pdls);

use threads;
use MCE::Shared;
use Time::HiRes 'time';

srand( 123 );

my $time = time;

my $n = 30000;      # input sample size
my $m = 10000;      # number of bootstrap repeats
my $r = $n;         # re-sample size

my $x   = random( $n );   $x->share_as('x');
my $avg = zeroes( $m );   $avg->share_as('avg');

my $seq = MCE::Shared->sequence( 0, $m - 1 );

sub parallel_task {
   srand;
   my ( $x, $avg ) = retrieve_pdls('x', 'avg');
   while ( defined ( my $seq_n = $seq->next() ) ) {
      my $idx  = random $r;
      $idx    *= $n;
      $avg->set( $seq_n, $x->index( $idx )->avg );
   }
}

threads->create( \&parallel_task ) for 1 .. 4;

# ... do other stuff ...

$_->join() for threads->list();

say $avg->pctover( pdl 0.05, 0.95 );
say time - $time, ' seconds';

__END__

# Output

[0.49395242 0.49936752]
1.28744792938232 seconds
[download]

Afterwards, re-validated PDL with MCE and released 1.847. The effort is mainly for folks running Perl lacking threads support. Here it is, PDL and MCE::Shared running similarly.

# https://www.perlmonks.org/?node_id=1214227

use strict;
use warnings;
use feature 'say';

use PDL;  # must load PDL before MCE::Shared

use MCE::Hobo;
use MCE::Shared 1.847;
use Time::HiRes 'time';

srand( 123 );

my $time = time;

my $n = 30000;      # input sample size
my $m = 10000;      # number of bootstrap repeats
my $r = $n;         # re-sample size

# On Windows, the non-shared piddle ($x) is unblessed in threads.
# Therefore, constructing the piddle inside the worker.
# UNIX platforms benefit from copy-on-write. Thus, one copy.

my $x   = ( $^O eq 'MSWin32' ) ? undef : random( $n );
my $avg = MCE::Shared->pdl_zeroes( $m );
my $seq = MCE::Shared->sequence( 0, $m - 1 );

sub parallel_task {
   $x = random( $n ) unless ( defined $x );

   while ( defined ( my $seq_n = $seq->next() ) ) {
      my $idx  = random $r;
      $idx    *= $n;
      # $avg is a shared piddle which resides inside the shared-
      # manager process or thread. The piddle is accessible via the
      # OO interface only.
      $avg->set( $seq_n, $x->index( $idx )->avg );
   }
}

MCE::Hobo->create( \&parallel_task ) for 1 .. 4;

# ... do other stuff ...

MCE::Hobo->wait_all();

# MCE sets the seed of the base generator uniquely between workers.
# Unfortunately, it requires running with one worker for predictable
# results (i.e. no guarantee in the order which worker computes the
# next input chunk).

say $avg->pctover( pdl 0.05, 0.95 );
say time - $time, ' seconds';

__END__

# Output

[0.49387191 0.49937053]
1.29038286209106 seconds
[download]

Regards, Mario

Here is the same thing using MCE. Workers obtain the next sequence number without involving the manager process. Thus, the reason why it runs faster. I had to think about it when I saw the run time.

# https://www.perlmonks.org/?node_id=1214227

use strict;
use warnings;
use feature 'say';

use PDL;  # must load PDL before MCE::Shared

use MCE 1.847;
use MCE::Shared 1.847;
use Time::HiRes 'time';

srand( 123 );

my $time = time;

my $n = 30000;      # input sample size
my $m = 10000;      # number of bootstrap repeats
my $r = $n;         # re-sample size

# On Windows, the non-shared piddle ($x) is unblessed in threads.
# Therefore, constructing the piddle inside the worker.
# UNIX platforms benefit from copy-on-write. Thus, one copy.

my $x   = ( $^O eq 'MSWin32' ) ? undef : random( $n );
my $avg = MCE::Shared->pdl_zeroes( $m );

MCE->new(
   max_workers => 4,
   sequence    => [ 0, $m - 1 ],
   chunk_size  => 1,
   user_begin  => sub {
      $x = random( $n ) unless ( defined $x );
   },
   user_func   => sub {
      my $idx  = random $r;
      $idx    *= $n;
      # $avg is a shared piddle which resides inside the shared-
      # manager process or thread. The piddle is accessible via the
      # OO interface only.
      $avg->set( $_, $x->index( $idx )->avg );
   }
)->run;

# MCE sets the seed of the base generator uniquely between workers.
# Unfortunately, it requires running with one worker for predictable
# results (i.e. no guarantee in the order which worker computes the
# next input chunk).

say $avg->pctover( pdl 0.05, 0.95 );
say time - $time, ' seconds';

__END__

# Output

[0.49387106  0.4993768]
1.09556317329407 seconds
[download]

Thank you, vr. I had no idea that PDL random is not unique between threads. MCE already sets the seed of the base generator, but did not do so for workers spawned as threads. This is resolved in MCE 1.847.

Regards, Mario

A reply from the future! Since 2.062, PDL's random numbers (using the random function) can be generated using pthreads since the new srand creates one seed per CPU.

This 'never to explicitly loop in vectorized language' answer, unfortunately, hides the ugly truth that for very large data we can end with huge R x M matrices of random indices and equally huge (equally unnecessary) matrices of all re-samplings, and thus die because of 'Out of memory!'.

Something else that would help here, as also discussed on #349, is more generalised first-class "index operations". Ideas and contributions here, on the GitHub issue, on the PDL mailing lists, or any other means are most welcome!

One other thought is that, for larger ndarrays (because POSIX threads have a startup cost), the use of vectorised operations is the way to harness multiple cores for free (for operations that support this), which a Perl for-loop cannot achieve.