PWC 248
PWC 248
I use PDL with Inline:Pdlpp for both challenges this week. I found these docs helpful: the perldocs on PDL::PP and PDL::API, the chapter on PP in the PDL book, and the Practical Magic guide to using PP.
Challenge 2 (submatrix sum)
We are given a matrix, say A, with M columns and N rows. We are asked to find a matrix, say B, with M-1 columns and N-1 rows, such that the entry for column i and row j of B, or B(i,j) using column-first notation is given by:
B(i,j) = A(i,j)+A(i,j+1)+A(i+1,j+1)+A(i+1,j)
With PP, we write XS code which uses a traditional pair of nested loops to create B. Here is the PP definition:
Pars => 'a(m,n); [o]b(p,q);',
RedoDimsCode => '$SIZE(p)=$SIZE(m)-1; $SIZE(q)=$SIZE(n)-1;',
Code => q{
PDL_Indx i, j, m_size, n_size;
m_size=$SIZE(m)-1;
n_size=$SIZE(n)-1;
for (i=0; i < m_size; i++) {
for (j=0; j < n_size; j++) {
$b(p=>i, q=>j)=
$a(m=>i, n=>j)+$a(m=>i+1, n=>j+1)+
$a(m=>i, n=>j+1)+$a(m=>i+1, n=>j);
}
}
},
);
The Pars section specifies the function signature. The submatrix_sum PDL method operates on a 2D PDL labelled a, with the two dimensions labelled m and n respectively. The output is also a 2D PDL, labelled b, with the two dimensions labelled p and q respectively.
We need to tell the compiler the size of B along dimensions p and q. We do this in the RedoDimsCode section.
Finally comes the main Code section. This is pretty much vanilla C with some additional custom preprocessor directives.
Challenge 1 (Shortest distance)
We are given a string say str and a single character say chr. We first find all the positions in the string where the character is chr. We then output an array showing for every position in the string, how far it is from the nearest occurrence of chr. For example, if the string is 'abaaaba' and the character is 'b', the output would be: [1,0,1,2,1,0,1]. The intermediate array showing the positions matching 'b' is [1,5]. (This is Perl, so arrays are zero-indexed).
My original plan was to set up two pp_def entries, one for an inner loop function that computes as output the sub-array showing the nearest elements, given as input a pair of consecutive items matching chr. For example, following the example above, the input could be the pair [1,5], and the output the sub-array 0,1,2,1,0 at positions 1 through 5 of the final output array.
I wanted to call this from another outer-loop function also defined in PP which calls the inner-loop function. This turned out to be very tricky and I was not able to get it done in time. I ended up just writing an inner-loop function and calling it from a PDL (Perl 5) outer loop, which is not very efficient. I also however wrote a single PP function that includes both the outer and inner loop, in a different script labelled ch-1a.pl.
Here is the pp_def for the original standalone inner-loop function.
'between',
Pars => 'indx a1(); indx a2(); int [o]b(n);',
RedoDimsCode => '$SIZE(n)=$a2()-$a1()+1;',
Code => q{
PDL_Indx i,c_mid, n_size=$SIZE(n);
if (0==(n_size % 2)){
c_mid=n_size/2;
}
else {
c_mid=(n_size+1)/2;
}
for (i=0; i < c_mid; i++) {
$b(n=>i)=i;
}
for (i=c_mid; i < n_size; i++) {
$b(n=>i)=n_size-i-1;
}
},
);
The logic reflects that if the gap between the two input indices leads to an odd-number-sized sub-array of distances, then the sub-array of distances would be of the pattern:
0,1,2,...x,x-1,...2,1,0
For example, for input-pair [1,5], we have [0,1,2,1,0].
If the gap leads to an even-number-sized sub-array of distances, then the sub-array of distances would be of the pattern:
0,1,2,...x-1,x,x,x-1,...2,1,0
For example, for input-pair [1,6], we have [0,1,2,2,1,0].
The XS function 'between' is called in PDL as a method, for example: pdl(1)->between(6).
I call it from a routine written in PDL (Perl 5):
my ($str,$chr)=@_;
($str =~ /$chr/) || return 0;
(1==length($chr)) || (die "Single character input needed:$!\n");
my $retval=zeros(length($str)); #-- pdl to hold return value
#-- get indices where the element equals chr
my @grp=(grep {$chr eq substr($str,$_,1)} 0..length($str));
#-- if the initial matching index is greater than 0 say x, fill 0:x
($grp[0] > 0) && ($retval->slice("0:$grp[0]") .= sequence($grp[0]+1)->slice('-1:0'));
#-- fill in distances between elements of @grp
map {
$retval->slice("$grp[$_]:$grp[$_+1]") .=
pdl($grp[$_])->between(pdl($grp[$_+1]));
} 0 .. ($#grp-1);
#-- if the last element of @grp is less than the length of $str, fill in the distance for the remaining elements.
($grp[-1] < (length($str)-1)) && ($retval->slice("$grp[-1]:") .= sequence(length($str)-$grp[-1]));
return $retval;
};
The array @grp contains the positions in str that match chr.
We loop through elements of @grp, calling the between operator at each stage. We also directly construct in Perl 5 the distances before the first index in @grp, as well as after the last index in @grp.
Potential improvements
If this was production code, I would tinker with it further, but given our one-week constraint, I will confine myself to reflecting on improvements in this blog.
Routines defined in PP should be capable of broadcasting to higher dimensions. My 'between' routine does mechanically broadcast, but not in a useful way. A better design could be to return a vector of the same length as the original string, with the distances inserted for the section between the input indices. For example, for input indices [1,6] with string length 10, the output would be [0,0,1,2,2,1,0,0,0,0,0], not just [0,1,2,2,1,0]. We could then say pass the routine more than one set of indices, say [[1,6],[6,8]] with the same string length, say 10, to get an output like
A more efficient approach is the one in ch-1a.pl, where I use a single pp_def function to take a list of indices and a string length as inputs, and return the shortest-distance array as output. Here is that single pp_def.
'shortest_distance_pp',
Pars => 'indx a(m); indx lenstr(); [o]b(n);',
RedoDimsCode => '$SIZE(n)=$lenstr();',
GenericTypes => ['N'],
Code => q{
PDL_Indx i,j,size_m, size_n, ctr_stop;
size_m=$SIZE(m);
size_n=$SIZE(n);
ctr_stop=$a(m=>0);
if (0 < $a(m=>0)){
for (i=0; i < ctr_stop; i++) {
$b(n=>i)=ctr_stop-i;
}
}
PDL_Indx a1,a2,c_mid,gap_size;
for (j=0; j<size_m; j++) {
a1=$a(m=>j);
a2=$a(m=>j+1);
gap_size=a2-a1+1;
if (0==(gap_size % 2)){
c_mid=gap_size/2;
}
else {
c_mid=(gap_size+1)/2;
}
for (i=a1; i < a1+c_mid; i++) {
$b(n=>i)=i-a1;
}
for (i=a1+c_mid; i <= a2; i++) {
$b(n=>i)=a2-i;
}
}
if (size_n > $a(m=>-1)+1){
a1=$a(m=>size_m-1);
for (i=a1; i < size_n; i++) {
$b(n=>i)=i-a1;
}
}
},
);
It follows the same logic as the Perl 5 outer loop + pp_def inner loop in ch-1.pl. But it would be faster and more scalable. Confining the input types to 'N' via the GenericTypes field ('N' refers to the PDL index type, usually a signed integer) reduces the work for the compiler which would otherwise create executable code unnecessarily for things like doubles etc.
Comments