Subtracting arrays of floats.

[Jonas Nilsson]

Hello.
I have an optimization problem.

I have two arrays of floats (millions of items in each). I need the absolut
sum of each pairwise difference.
I would prefer to have the arrays as packed arrays (double or single
precision will both work) since I would otherwise need to much memory.

I need to do this quite quickly. The code that i now use takes about 5 days
to run on a 2GHz pentium. Is there any trick, with bitwise operators or some
other tweak what would significantly speed op this code? I'm sad so say that
I don't have the knowledge about using C.

/jN

Some example code just for you to see what I wan't to acheve.
use strict;
use Benchmark;

my $count=1e5; #Just for testing

my @array1=map {rand()-0.5} (1..$count);
my @array2=map {rand()-0.5} (1..$count);
my $array1_f=pack "f*",@array1;
my $array2_f=pack "f*",@array2;
my $array1_d=pack "d*",@array1;
my $array2_d=pack "d*",@array2;


timethese(20,{
'array'=>sub {
my $abssum=0;
for (0..$count-1) {
$abssum+=abs($array1[$_]-$array2[$_]);
}
},
'scalar of f'=>sub {
my $abssum=0;
for (0..$count-1) {

$abssum+=abs(unpack('f',substr($array1_f,$_*4,4))-unpack('f',substr($array2_
f,$_*4,4)));
}
},
'scalar of d'=>sub {
my $abssum=0;
for (0..$count-1) {

$abssum+=abs(unpack('d',substr($array1_d,$_*8,8))-unpack('d',substr($array2_
d,$_*8,8)));
}
},

});

 

[Tassilo v. Parseval]

Also sprach Jonas Nilsson:

I have an optimization problem.

Rather a language problem.

I have two arrays of floats (millions of items in each). I need the absolut
sum of each pairwise difference.
I would prefer to have the arrays as packed arrays (double or single
precision will both work) since I would otherwise need to much memory.

I need to do this quite quickly. The code that i now use takes about 5 days
to run on a 2GHz pentium. Is there any trick, with bitwise operators or some
other tweak what would significantly speed op this code? I'm sad so say that
I don't have the knowledge about using C.

This is the right moment to learn it. Honestly, the real problem here is
the choice of tools. Perl is not the right one. Here's a solution
written in XS (that is the C dialect used to create Perl modules).
Create the module stub like this:

h2xs -A -b 5.5.3 -n Math:airwiseDiff

Then delve into Math/PairwiseDiff and modify PairwiseDiff.xs. Here's a
solution that should be performant enough for a start. Your .xs file
should look like this:

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include <math.h>

#include "ppport.h"

#define ptr2ary(a) (ARY*) SvIV((SV*)SvRV(a))

typedef struct ARY {
float *f;
int num;
} ARY;

MODULE = Math:airwiseDiff PACKAGE = Math:airwiseDiff

void
new (CLASS, ...)
char *CLASS;
PREINIT:
ARY *ary;
register float *s;
register int i;
SV *ret;
PPCODE:
New(0, ary, 1, ARY);
New(0, ary->f, items - 1, float);
s = ary->f;
for (i = 1; i < items; i++)
*s++ = SvNV(ST(i));
ary->num = items;
ret = sv_newmortal();
/* bless a pointer to the C struct into this class */
sv_setref_pv(ret, CLASS, (void*)ary);
XPUSHs(ret);
XSRETURN(1);

void
compute_diff (ary1, ary2)
SV *ary1;
SV *ary2;
PREINIT:
float diff;
register ARY *a1, *a2;
int num;
register int i;
SV *ret;
float *s1, *s2; /* start-pointers of floats */
PPCODE:
/* turn Perl objects into C structs */
a1 = ptr2ary(ary1);
a2 = ptr2ary(ary2);

diff = 0;
num = a1->num;
s1 = a1->f;
s2 = a2->f;
for (i = 0; i < num; i++)
diff += fabsf(*s1++ - *s2++);

ret = newSVnv(diff);
XPUSHs(sv_2mortal(ret));
XSRETURN(1);

void
DESTROY (ary)
SV *ary;
PREINIT:
ARY *a;
PPCODE:
a = ptr2ary(ary);
Safefree(a->f);
Safefree(a);

I omitted some otherwise maybe necessary checks and typecasts. 'float'
could be replaced with 'unsigned float' in some cases, or even 'double'.
Should work nonetheless as long as the two arrays have the same number
of elements (it's not too hard to make it work for other cases as well).

The Makefile.PL needs to have an additional linker target:

...
'LIBS' => ['-lm'],
...

The Perl code using the above (after 'perl Makefile.PL; make etc.' looks
like this:

use Math:airwiseDiff;

my $a1 = Math:airwiseDiff->new( map { rand - 0.5 } 0 .. 200_000 );
my $a2 = Math:airwiseDiff->new( map { rand - 0.5 } 0 .. 200_000 );

print $a1->compute_diff($a2), "\n";

# just for testing: should be the same value
print $a2->compute_diff($a1);

The slowest part is the constructor new() because it has to convert Perl
scalar variables into C floats. It takes 1 or 2 secs on my machine
(Athlon 900) for the above 200_001 floats.

compute_diff() is then very fast. For a million plus 1 floats in each
array, the runtime of the whole script is:

ethan@ethan:~/Projects/Math/PairwiseDiff$ time perl test.pl
777.050897435469
777.050897435469
real 0m4.784s
user 0m4.420s
sys 0m0.110s

So it's less than 5 seconds. It'll be faster on your computer. It
shouldn't eat too much memory either.

There are some Perl modules on the CPAN that you could use instead. PDL
could likewise be quite fast since it's written in C as well. Doing it
on your own however has the advantage that you can make some
optimizations matching exactly your given problem.

Tassilo

 

[Anno Siegel]

Hello.
I have an optimization problem.

I have two arrays of floats (millions of items in each). I need the absolut
sum of each pairwise difference.

That looks like a problem for the Perl Data Language, available as
PDL from CPAN. I can't offer an immediate solution -- I never got
around to using PDL myself -- but that's the kind of problem it was
designed for.

An individual XS solution, as Tassilo suggests in another reply may
be the absolute fastest one attainable. If you use the Inline module,
it doesn't have to look quite as formidable. See my followup to
Tassilo's post.

[useful explanation of problem snipped]

Anno

 

[Anno Siegel]

Also sprach Jonas Nilsson:


Rather a language problem.


This is the right moment to learn it. Honestly, the real problem here is
the choice of tools. Perl is not the right one. Here's a solution
written in XS (that is the C dialect used to create Perl modules).
Create the module stub like this:

h2xs -A -b 5.5.3 -n Math:airwiseDiff

[good instructions how to build the XS module snipped]

The Inline module makes simple things simple. It knows a lot about
XS that the user doesn't have to know (at least in simple cases like
this one). Here is how a similar function (named "distance" here) can
be set up using Inline. We pass in two strings that contain the packed
doubles, plus the number of data points. The latter is just for
convenience, we could determine the length of one of the strings and
use that.


#!/usr/bin/perl
use strict; use warnings; $| = 1; # @^~`

# number of data points
use constant N => 100_000;

# set up the data arrays
my ( $array1_d, $array2_d);
for ( 0 .. N-1 ) {
$array1_d .= pack 'd', rand() - 0.5;
$array2_d .= pack 'd', rand() - 0.5;
}

my $res = distance( $array1_d, $array2_d, N);
print "$res\n";

use Inline C => <<EOC;

double distance( char* x, char* y, int n) {
/* cast pointers */
double* xx = (double *) x;
double* yy = (double *) y;
/* build sum */
double sum = 0.0;
while ( n-- ) {
sum += fabs( *xx++ - *yy++);
}
return( sum);
}
EOC
__END__

Anno

 

[Jonas Nilsson]

The Inline module makes simple things simple. It knows a lot about
XS that the user doesn't have to know (at least in simple cases like
this one). Here is how a similar function (named "distance" here) can
be set up using Inline. We pass in two strings that contain the packed
doubles, plus the number of data points. The latter is just for
convenience, we could determine the length of one of the strings and
use that.

Ok. Call me thick, but I'm not an expert, to say the least.

To compile C either in Inline module or as a separate XS-module i guess I
need a C-compiler. So I don't have one... The nmake utility for Windows
isn't a C-compiler or is it?

When I try to run/make code in one of the above ways I get DOS complaining
about the lack of cl command. Is this a compiler?

What exactly do I need to run this?
I've got:
1. A computer
2. An operating system (Windows 2000)
3. perl5.6.1 from activestate
4. nmake15 (from windows support) for WinNT/Win95. I couldn't find one for
Win2000.
5. The Inline module via ppm from
http://ppm.ActiveState.com/cgibin/PPM/ppmserver.pl?urn:/PPMServer

And nothing more (...well I've got a nice espresso machine, but I doubt it
would influence much).

When I try to use the Inline C i get from DOS:
'cl' is not recognized as an internal or external command, operable program
or batch file. (and a lot more)

Please tell me what I need.
/jN

 

[Tassilo v. Parseval]

Also sprach Anno Siegel:

This is the right moment to learn it. Honestly, the real problem here is
the choice of tools. Perl is not the right one. Here's a solution
written in XS (that is the C dialect used to create Perl modules).
Create the module stub like this:

h2xs -A -b 5.5.3 -n Math:airwiseDiff

[good instructions how to build the XS module snipped]

The Inline module makes simple things simple. It knows a lot about
XS that the user doesn't have to know (at least in simple cases like
this one). Here is how a similar function (named "distance" here) can
be set up using Inline. We pass in two strings that contain the packed
doubles, plus the number of data points. The latter is just for
convenience, we could determine the length of one of the strings and
use that.

Very terse, nice! Nonetheless I keep thinking that Inline::C only makes
those things simple that are simple by nature (and therefore simple in
XS).

But it looks less scary which is a good thing for someone starting with C
extensions for Perl.

#!/usr/bin/perl
use strict; use warnings; $| = 1; # @^~`

# number of data points
use constant N => 100_000;

# set up the data arrays
my ( $array1_d, $array2_d);
for ( 0 .. N-1 ) {
$array1_d .= pack 'd', rand() - 0.5;
$array2_d .= pack 'd', rand() - 0.5;
}

my $res = distance( $array1_d, $array2_d, N);
print "$res\n";

use Inline C => <<EOC;

double distance( char* x, char* y, int n) {
/* cast pointers */
double* xx = (double *) x;
double* yy = (double *) y;

Using pack() and these casts is quite smart. Besides being a fast way
it'll also use very little memory.

In my solution, 1_000_000 floats stored in 'ARY*' need around 480kB
memory (at least on machines with a memory-alignment of 4 or 8 bytes
which should be common enough). So that's ok. But perl still has to
build a million scalars and put them onto the stack. And since perl
will only reluctantly give the memory back to the OS (if ever), this
could be a real problem.

/* build sum */
double sum = 0.0;
while ( n-- ) {
sum += fabs( *xx++ - *yy++);
}
return( sum);
}
EOC
__END__

One minor thing: you use fabs() from math.h. Is Inline::C smart enough
to automatically link against libm? It might be necessary to add a

use Inline C => Config =>
LIBS => '-lm';

Tassilo

 

[Tassilo v. Parseval]

Also sprach Jonas Nilsson:

Ok. Call me thick, but I'm not an expert, to say the least.

To compile C either in Inline module or as a separate XS-module i guess I
need a C-compiler. So I don't have one... The nmake utility for Windows
isn't a C-compiler or is it?

Not having a C compiler is a little problem now. I heard rumours that
there is a free version of VisualC available but I couldn't find it.
nmake is not a compiler, it's just the make tool from Microsoft (which
is needed as well).

When I try to run/make code in one of the above ways I get DOS complaining
about the lack of cl command. Is this a compiler?

Yes, this is the compiler you need.

What exactly do I need to run this?
I've got:
1. A computer
2. An operating system (Windows 2000)
3. perl5.6.1 from activestate
4. nmake15 (from windows support) for WinNT/Win95. I couldn't find one for
Win2000.
5. The Inline module via ppm from
http://ppm.ActiveState.com/cgibin/PPM/ppmserver.pl?urn:/PPMServer

And nothing more (...well I've got a nice espresso machine, but I doubt it
would influence much).

Quite.

The thing that you are lacking is indeed the compiler. Since the code to
solve your problem already exists in this group, someone would just have
to make a PPM distribution out of it. Maybe I'll boot Windows later
today and can do that (either Anno's or my version, but probably a
healthy mix of both...I liked the pack() so much).

Hmmh, ActiveState's new automatic build-process could be abused for
that, I guess. Create a Perl/XS/Inline module and upload it to the CPAN.
ActiveState currently picks up every CPAN upload and tries to convert it
into a PPM distribution. This includes modules in need of a C compiler.
Once it has made it into the official PPM repositories, the module could
be deleted again from the CPAN.

With a little bit of Makefile.PL trickery, one could possibly make them
compile Win32 applications for oneself.

Tassilo

 

[Anno Siegel]

Also sprach Anno Siegel:

Tassilo v. Parseval <[email protected]> wrote in comp.lang.perl.misc:

[...]

/* build sum */
double sum = 0.0;
while ( n-- ) {
sum += fabs( *xx++ - *yy++);
}
return( sum);
}
EOC
__END__

One minor thing: you use fabs() from math.h. Is Inline::C smart enough
to automatically link against libm? It might be necessary to add a

use Inline C => Config =>
LIBS => '-lm';

I don't know, maybe I just got lucky, but it found fabs() for me. It may
well be that the library must be given in a different setup.

Anno

 

[Tassilo v. Parseval]

Also sprach Anno Siegel:

Also sprach Anno Siegel:

Tassilo v. Parseval <[email protected]> wrote in comp.lang.perl.misc:

[...]

/* build sum */
double sum = 0.0;
while ( n-- ) {
sum += fabs( *xx++ - *yy++);
}
return( sum);
}
EOC
__END__

One minor thing: you use fabs() from math.h. Is Inline::C smart enough
to automatically link against libm? It might be necessary to add a

use Inline C => Config =>
LIBS => '-lm';

I don't know, maybe I just got lucky, but it found fabs() for me. It may
well be that the library must be given in a different setup.

Hmmh, right. My XS thing also works without it. I thought that maybe the
linker flags in Config.pm were copied but this doesn't happen:

[...]
cc -c -I. -O3 -march=athlon -fno-strict-aliasing -D_LARGEFILE_SOURCE
-D_FILE_OFFSET_BITS=64 -O2 -DVERSION=\"0.01\" -DXS_VERSION=\"0.01\"
-fpic "-I/usr/lib/perl5/5.8.0/i686-linux/CORE" PairwiseDiff.c
Running Mkbootstrap for Math:airwiseDiff ()
chmod 644 PairwiseDiff.bs
rm -f blib/arch/auto/Math/PairwiseDiff/PairwiseDiff.so
LD_RUN_PATH="" cc -shared -L/usr/local/lib PairwiseDiff.o -o
blib/arch/auto/Math/PairwiseDiff/PairwiseDiff.so

Anyway, one shouldn't complain when things unexpectedly work.

Tassilo

 

[Tassilo v. Parseval]

[ posted and mailed ]

Also sprach Jonas Nilsson:

When I try to use the Inline C i get from DOS:
'cl' is not recognized as an internal or external command, operable program
or batch file. (and a lot more)

You are lucky. I just incorporated Anno's pack() trick into my module
and it compiled without any problems under Windows/VisualC (this is
rather the exception, not the rule). I made a .ppd for it which you can
now install even if no VisualC is around [ might line-wrap ]:

C:\> ppm install http://www-users.rwth-aachen.de/tassilo.parseval/Math-PairwiseDiff.ppd

I don't have WinNT/2000 but only WinME. It should nonetheless work for
you.

It comes with some rudimentary POD-documentation as well. This little
script demonstrates all you can do with it (and yes, it still ran under
Windows):

#! /usr/bin/perl

use strict;
use Math:airwiseDiff;
use Data:umper;

my ($ary1_d, $ary2_d);

$ary1_d .= pack('d', rand() - 0.5) for 0 .. 100 ;
$ary2_d .= pack('d', rand() - 0.5) for 0 .. 100 ;

my $a1 = Math:airwiseDiff->new($ary1_d);
my $a2 = Math:airwiseDiff->new($ary2_d);

print Dumper [ $a1->to_ary ];
print "Num of elements: ", scalar $a1->to_ary, "\n";

print $a1->diff($a2), "\n";
print $a2->diff($a1);

So I just added the to_ary() method that returns the packed doubles as a
list of real Perl scalars containing human-readable floating point
numbers. It now also works with arrays of different length.
Furthermore, you no longer have to pass the number of doubles in the
packed string to the constructor. It can figure that out itself.

This module is not generic enough to be added to the CPAN however.

Tassilo

 

[Tassilo v. Parseval]

Also sprach Tassilo v. Parseval:

[ posted and mailed ]

Or rather, posted and tried to mail. Naturally, the email address
bounced.

Also sprach Jonas Nilsson:

Tassilo

 

[Jonas Nilsson]

[ posted and mailed ]

Also sprach Jonas Nilsson:

When I try to use the Inline C i get from DOS:
'cl' is not recognized as an internal or external command, operable program
or batch file. (and a lot more)

You are lucky. I just incorporated Anno's pack() trick into my module
and it compiled without any problems under Windows/VisualC (this is
rather the exception, not the rule). I made a .ppd for it which you can
now install even if no VisualC is around [ might line-wrap ]:

Well, thank you a million times. Sometimes it's hard to believe that there
really are people that would help a complete stranger for free...

Could you please post the revised source also, for my inspection. I would
like to try making some modifications. I'll se if I can find any free
compiler. Maybee I'll start learing C.

If I were to get a complier, what would my specifications be. Are they all
the same? I mean, I just need the complier, not a whole bundle of fancy
editors and such.

Thanks again.
/jN

 

[Tassilo v. Parseval]

Also sprach Jonas Nilsson:

[ posted and mailed ]

Also sprach Jonas Nilsson:

When I try to use the Inline C i get from DOS:
'cl' is not recognized as an internal or external command, operable
program or batch file. (and a lot more)

You are lucky. I just incorporated Anno's pack() trick into my module
and it compiled without any problems under Windows/VisualC (this is
rather the exception, not the rule). I made a .ppd for it which you can
now install even if no VisualC is around [ might line-wrap ]:

Well, thank you a million times. Sometimes it's hard to believe that there
really are people that would help a complete stranger for free...

Could you please post the revised source also, for my inspection. I would
like to try making some modifications. I'll se if I can find any free
compiler. Maybee I'll start learing C.

Good idea. The (shortened) link to the source is

<http://xrl.us/th5>

Yet, in order to understand and modify it just knowing C might not be
enough. XS code uses the perlapi with a lot of funny looking macros and
such. This is all documented:

perlguts.pod
perlxstut.pod
perlxs.pod
perlapi.pod
perlclib.pod # describes the perlapi versions of libc functions

I used these symbols from the perlapi (in order of appearance):

SvPV
New
Copy
sv_newmortal
sv_setref_pv
XPUSHs
XSRETURN
SvIV
SvRV
newSVnv
GIMME
G_SCALAR
sv_2mortal
G_ARRAY
EXTEND
PUSHs
Safefree

either to be found in perlapi or perlclib (New, Copy, Safefree). So you
can look them up while trying to make sense of the code.

I learned C through XS, actually. You can do it with Inline::C if you
prefer. The above manpages are then still valuable and worth a read.

If I were to get a complier, what would my specifications be. Are they all
the same? I mean, I just need the complier, not a whole bundle of fancy
editors and such.

Usually just a compiler-suit. The compiler itself has other necessary
components, the linker for instance. This is usually included since
otherwise the compiler wouldn't be of much use. And a make tool. nmake
will be ok if you use VisualC.

Since you are doing this under Windows, maybe you have a look at

<http://www.cygwin.com>

It's a Linux-environment running under Windows. You have a proper shell
and all the essential tools to develop software. The advantage is that
it all fits together very well (and is free, of course). You get a real
make program plus the gcc as C compiler. cygwin has software packages
that you can selectively install. So you'd probably want to install
perl, make, gcc, maybe a good text-editor (although you can still use
your current one) and some related libraries.

I think with a little bit of tweaking, you can even use this cygwin
environment to compile modules for ActivePerl. I remember that I tried
it a couple of times and it usually worked pretty well.

Or you just install Linux, which is the most convenient way.

Anyway, if you can get hold of VisualC from whatever legal or illegal
sources (not my concern, you don't need cygwin. In the beginning it's
a little harder since Windows does things differently from cygwin (which
does it the UNIX-way). Also, perlxstut (the XS tutorial) is assuming a
UNIXish environment. So it's easier to follow and experiment with the
examples.

Tassilo