Processing values of a piddle (PDL) speedup using 'at' vs. 'index'

[kevbot]

Hello Monks,

I have been using PDL in some data processing scripts that I use at work. One of my modules takes float values from a pdl and processes them before printing them to a text file. Although a text file is very inefficient for the type of data I am handling, I need to output a text file for use with some legacy (proprietary) software.

While profiling the code that writes the text file (thank you Devel::NYTProf), I discovered that there was a lot of time spent in various "PDL" subroutines. It occurred to me that I had made a bad assumption (or perhaps a case of expecting PDL to DWIM but not matching my expectations). I thought this issue might be encountered by others; hence, I'm writing this node.

I have been happy to see that PDL has been getting some love recently (new version released in May 2012). Also, Joel Berger has been making some blog posts recently regarding using Perl and PDL in scientific applications.

OK, back to the problem at hand. I had assumed (incorrectly) that when I placed a single value from a 1-dimensional piddle into a variable that the resulting value would be a perl scalar. For example,

my $val = $my_pdl->index($i);

In reality the $val is a piddle (i.e. ref($pdl) returns 'PDL'). My code still ran and gave the expected results...however it ran more slowly (compared to the other method using at, see below).

In my use case, I take a pre-existing pdl that contains hundreds of thousands (perhaps millions) of float values and check each element one-by-one (e.g. to see if it equals a special value, etc.) and then change the resulting value before I put it back into a perl array...where I latter dump the array into a text file. It turns out that the code runs significantly faster if I simply change the code to,

my $val = $my_pdl->at($i);

According to the PDL docs, the at method returns a single value inside a piddle as perl scalar (whereas, the index method returns a pdl). Perhaps this is a long write-up for a simple case of confusion, but the speed difference was quite noticeable. The code below gives a benchmark of iterating over all elements of a 10,000 element pdl using both methods (and doing nothing else).

#!/usr/bin/env perl

use strict;
use warnings;

use PDL;
use Benchmark qw(:all);

cmpthese(
    100,
    {   'pdl_values'  => sub {&pdl_value},
        'perl_values' => sub {&perl_scalar_value},
    }
);

exit;

sub pdl_value {
    my $pdl = ones( float, 10000 );
    my $nelem = nelem($pdl);

    for ( my $i = 0; $i < $nelem; ++$i ) {
        my $val = $pdl->index($i);

        #Do something
    }

    return;
}

sub perl_scalar_value {
    my $pdl           = ones( float, 10000 );
    my $nelem         = nelem($pdl);

    for ( my $i = 0; $i < $nelem; ++$i ) {
        my $val = $pdl->at($i);

        #Do something
    }

    return;
}

              Rate  pdl_values perl_values
pdl_values  16.8/s          --        -74%
perl_values 64.1/s        282%          --

When I add the code to put the processed values into a perl array and to check for a "special value" the difference is even larger.

#!/usr/bin/env perl

use strict;
use warnings;

use PDL;
use Benchmark qw(:all);

cmpthese(
    100,
    {   'pdl_values'  => sub {&pdl_value},
        'perl_values' => sub {&perl_scalar_value},
    }
);


exit;

sub pdl_value {
    my $pdl = ones( float, 10000 );
    $pdl->index(0) .= 999;
    my $nelem         = nelem($pdl);
    my $special_value = 999;

    my @values;
    for ( my $i = 0; $i < $nelem; ++$i ) {
        my $val = $pdl->index($i);
        if ( $val == $special_value ) {
            $val = undef;
        }
        push @values, $val;
    }

    #Do something with @values
    return;
}

sub perl_scalar_value {
    my $pdl = ones( float, 10000 );
    $pdl->index(0) .= 999;
    my $nelem         = nelem($pdl);
    my $special_value = 999;

    my @values;
    for ( my $i = 0; $i < $nelem; ++$i ) {
        my $val = $pdl->at($i);
        if ( $val == $special_value ) {
            $val = undef;
        }
        push @values, $val;
    }

    #Do something with @values
    return;
}

              Rate  pdl_values perl_values
pdl_values  1.51/s          --        -97%
perl_values 52.1/s       3347%          --

I discovered the at method by reading the PDL Book. I highly recommend this book for those that are just getting started with PDL. It can be a little difficult to find what you are looking for in the PDL docs (compared to the book). The book actually warns against the use of at, since it is "slow"; however, it is a significant improvement in this case. Perhaps there is yet another way that would be faster. Here's the relevant quote from the PDL Book:

Conversion to Perl types: at and list

You can get a PDL scalar out into the Perl world with at, which requires the index of the scalar to pull out:

pdl> $a = xvals(5)*2; # $a is a PDL
pdl> $a4 = $a->at(4); # $a4 is a perl scalar

You can also export a whole PDL with list:

pdl> @a = $a->list; pdl> for($a->list) { print $_, - ; }
0-2-4-6-8-

Be careful with at, as you almost never want to use it - it is tedious for anything nontrivial, and extremely slow! Particularly if you find yourself placing an at call inside a for loop, you should probably stop and think about how to use threading for your problem - see below.

Well, I'm off to learn a little more about PDL.

UPDATE: I removed an extraneous line of code from the first code example (it didn't affect the benchmark results).