Task 2: Submatrix Sum
Submitted by: Jorg Sommrey
You are given a NxM matrix A of integers.
Write a script to construct a (N-1)x(M-1) matrix B
having elements that are the sum over the 2x2 submatrices of A,
b[i,k] = a[i,k] + a[i,k+1] + a[i+1,k] + a[i+1,k+1]
Example 1
Input: $a = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12]
]
Output: $b = [
[14, 18, 22],
[30, 34, 38]
]
####
sub sms_WxH_Perl_4loops( $m, $w, $h ) {
my ( $W, $H ) = ( scalar $m-> [0]-> @*, scalar @$m );
my @ret;
for my $i ( 0 .. $H - $h ) {
for my $j ( 0 .. $W - $w ) {
for my $k ( 0 .. $w - 1 ) {
for my $l ( 0 .. $h - 1 ) {
$ret[ $i ][ $j ] += $m-> [ $i + $l ][ $j + $k ];
}
}
}
}
return \@ret
}
sub sms_WxH_Perl_3loops( $m, $w, $h ) {
my ( $W, $H ) = ( scalar $m-> [0]-> @*, scalar @$m );
my ( @tmp, @ret );
for my $i ( 0 .. $H - 1 ) {
for my $j ( 0 .. $W - $w ) {
for my $k ( 0 .. $w - 1 ) {
$tmp[ $i ][ $j ] += $m-> [ $i ][ $j + $k ];
}
next if $i < $h - 1;
for my $k ( 0 .. $h - 1 ) {
$ret[ $i - $h + 1 ][ $j ] += $tmp[ $i - $k ][ $j ];
}
}
}
return \@ret
}
$m = [
map [
map rand, 0 .. 149
], 0 .. 149
];
__END__
+
Time (s) vs. N (NxN submatrix, 150x150 matrix)
+----------------------------------------------------------+
|+ + + + + + + + |
10 |-+ +-|
|+ +|
|+ +|
|+ A A A A +|
|+ A A +|
|+ A A +|
| A A |
| |
1 |-+ A +-|
|+ A +|
|+ +|
|+ +|
|+ A +|
| A |
|+ +|
| B B B B B B |
0.1 |-+ B B B +-|
|+ B +|
|+ B B +|
|+ B +|
|+ +|
| |
|+ B +|
| |
0.01 |-+ +-|
|+ + + + + + + + +|
+----------------------------------------------------------+
0 20 40 60 80 100 120 140
sms_WxH_Perl_4loops A
sms_WxH_Perl_3loops B
+-----+-------+-------+
| N | A | B |
+-----+-------+-------+
| 10 | 0.250 | 0.053 |
| 20 | 0.809 | 0.091 |
| 30 | 1.528 | 0.116 |
| 40 | 2.288 | 0.134 |
| 50 | 2.988 | 0.144 |
| 60 | 3.512 | 0.147 |
| 70 | 3.828 | 0.144 |
| 80 | 3.862 | 0.134 |
| 90 | 3.628 | 0.122 |
| 100 | 3.122 | 0.103 |
| 110 | 2.453 | 0.088 |
| 120 | 1.675 | 0.066 |
| 130 | 0.900 | 0.044 |
| 140 | 0.284 | 0.022 |
+-----+-------+-------+
####
sub sms_2x2_PDL ( $m ) {
$m-> slice( '0:-2', '0:-2' ) +
$m-> slice( '1:-1', '0:-2' ) +
$m-> slice( '0:-2', '1:-1' ) +
$m-> slice( '1:-1', '1:-1' )
}
####
sub sms_WxH_PDL_naive ( $m, $w, $h ) {
my $ret = 0;
for my $r ( 0 .. $h - 1 ) {
my $r_sel = "$r:" . ( $r - $h );
for my $c ( 0 .. $w - 1 ) {
$ret += $m-> slice( "$c:" . ( $c - $w ), $r_sel )
}
}
return $ret
}
####
use Inline Pdlpp => <<'EOPP';
pp_def('sms_WxH_pdlpp_4loops',
Pars => 'a(m,n); [o]b(p,q);',
OtherPars => 'PDL_Indx par_sm_w; PDL_Indx par_sm_h;',
RedoDimsCode => '$SIZE(p) = $SIZE(m) - $COMP(par_sm_w) + 1;
$SIZE(q) = $SIZE(n) - $COMP(par_sm_h) + 1;',
Code => '
PDL_Indx i, j, k, l, w, h, W, H;
w = $COMP(par_sm_w);
h = $COMP(par_sm_h);
W = $SIZE(m);
H = $SIZE(n);
for (i = 0; i < H - h + 1; i++) {
for (j = 0; j < W - w + 1; j++) {
for (k = 0; k < w; k++) {
for (l = 0; l < h; l++) {
$b(p=>j, q=>i) += $a(m=>j+k, n=>i+l);
}
}
}
}
',
);
pp_def('sms_WxH_pdlpp_3loops',
Pars => 'a(m,n); [o]b(p,q); [t]t(p,n)',
OtherPars => 'PDL_Indx par_sm_w; PDL_Indx par_sm_h;',
RedoDimsCode => '$SIZE(p) = $SIZE(m) - $COMP(par_sm_w) + 1;
$SIZE(q) = $SIZE(n) - $COMP(par_sm_h) + 1;',
Code => '
PDL_Indx i, j, k, w, h, W, H;
w = $COMP(par_sm_w);
h = $COMP(par_sm_h);
W = $SIZE(m);
H = $SIZE(n);
for (i = 0; i < H; i++) {
for (j = 0; j < W - w + 1; j++) {
for (k = 0; k < w; k++) {
$t(p=>j, n=>i) += $a(m=>j+k, n=>i);
}
if (i >= h - 1) {
for (k = 0; k < h; k++) {
$b(p=>j, q=>i-h+1) += $t(p=>j, n=>i-k);
}
}
}
}
',
);
EOPP
####
sub sms_WxH_PDL_lags ( $m, $w, $h ) {
my ( $W, $H ) = $m-> dims;
return $m
-> lags( 0, 1, 1 + $W - $w )
-> sumover
-> lags( 1, 1, 1 + $H - $h )
-> xchg( 0, 1 )
-> sumover
-> slice( '-1:0,-1:0' )
}
$m = 1 + sequence 4, 3;
print $m,
sms_WxH_PDL_lags( $m, 2, 2 ),
sms_WxH_PDL_lags( $m, 3, 3 ),
sms_WxH_PDL_lags( $m, 4, 3 ),
sms_WxH_PDL_lags( $m, 1, 1 );
__END__
[
[ 1 2 3 4]
[ 5 6 7 8]
[ 9 10 11 12]
]
[
[14 18 22]
[30 34 38]
]
[
[54 63]
]
[
[78]
]
[
[ 1 2 3 4]
[ 5 6 7 8]
[ 9 10 11 12]
]
####
+
Time (s) vs. N (NxN submatrix, PDL: Double D [1500,1500] matrix)
+-----------------------------------------------------------+
2 |-+ + + + + + + + + +-|
| C C |
| C |
| C |
| D D D D C |
1.5 |-+ A C +-|
| D D |
| C |
| C D |
| |
| B D D |
1 |-+ +-|
| C |
| C |
| D D |
| |
| A D |
0.5 |-+ C +-|
| D |
| |
| B D |
| |
| DC |
0 |-+ DD +-|
| + + + + + + + + |
+-----------------------------------------------------------+
0 200 400 600 800 1000 1200 1400
sms_WxH_PDL_naive A
sms_WxH_pdlpp_4loops B
sms_WxH_pdlpp_3loops C
sms_WxH_PDL_lags D
+------+-------+-------+-------+-------+
| N | A | B | C | D |
+------+-------+-------+-------+-------+
| 2 | 0.038 | 0.009 | 0.009 | 0.047 |
| 4 | 0.294 | 0.031 | 0.019 | 0.041 |
| 6 | 0.572 | 0.072 | 0.019 | 0.044 |
| 10 | 1.531 | 0.250 | 0.031 | 0.031 |
| 20 | | 1.150 | 0.072 | 0.022 |
| 100 | | | 0.494 | 0.394 |
| 200 | | | 0.931 | 0.797 |
| 300 | | | 1.263 | 1.153 |
| 400 | | | 1.544 | 1.431 |
| 500 | | | 1.756 | 1.594 |
| 600 | | | 1.869 | 1.653 |
| 700 | | | 1.872 | 1.641 |
| 800 | | | 1.788 | 1.594 |
| 900 | | | 1.597 | 1.481 |
| 1000 | | | 1.362 | 1.288 |
| 1100 | | | 1.144 | 1.081 |
| 1200 | | | 0.841 | 0.809 |
| 1300 | | | 0.575 | 0.572 |
| 1400 | | | 0.291 | 0.272 |
+------+-------+-------+-------+-------+
####
sub _do_dimension ( $m, $w ) {
$m -> slice([ 0, $w - 1 ], [])
-> sumover
-> slice( '*' )
-> append(
$w >= ( $m-> dims )[ 0 ]
? empty
: $m-> slice([ $w, -1 ], []) -
$m-> slice([ 0, -1 - $w ], []))
-> cumusumover
-> xchg( 0, 1 )
}
sub sms_WxH_PDL_sliding ( $m, $w, $h ) {
_do_dimension
_do_dimension( $m, $w ), $h
}
$m = 1 + sequence 4, 3;
print $m,
sms_WxH_PDL_sliding( $m, 2, 2 ),
sms_WxH_PDL_sliding( $m, 3, 3 ),
sms_WxH_PDL_sliding( $m, 4, 3 ),
sms_WxH_PDL_sliding( $m, 1, 1 );
__END__
# same (correct) output as above, skipped
+
Time (s) vs. N (NxN submatrix, PDL: Double D [1500,1500] matrix)
+-----------------------------------------------------------+
| + + + + + + + + |
| A |
| A |
| A A |
1.5 |-+ A +-|
| A |
| |
| A |
| |
| A |
| A |
1 |-+ +-|
| |
| A |
| A |
| |
| |
| A |
0.5 |-+ +-|
| |
| A |
| A |
| |
| B B |
| BB B B B B B B B B B |
0 |-+ B B B B +-|
| + + + + + + + + |
+-----------------------------------------------------------+
0 200 400 600 800 1000 1200 1400
sms_WxH_PDL_lags A
sms_WxH_PDL_sliding B
+------+-------+-------+
| N | A | B |
+------+-------+-------+
| 2 | 0.025 | 0.100 |
| 4 | 0.034 | 0.075 |
| 6 | 0.016 | 0.084 |
| 10 | 0.034 | 0.091 |
| 20 | 0.028 | 0.084 |
| 100 | 0.369 | 0.081 |
| 200 | 0.800 | 0.078 |
| 300 | 1.122 | 0.116 |
| 400 | 1.416 | 0.050 |
| 500 | 1.572 | 0.053 |
| 600 | 1.644 | 0.022 |
| 700 | 1.675 | 0.025 |
| 800 | 1.587 | 0.019 |
| 900 | 1.456 | 0.019 |
| 1000 | 1.297 | 0.009 |
| 1100 | 1.047 | 0.009 |
| 1200 | 0.813 | 0.019 |
| 1300 | 0.566 | 0.006 |
| 1400 | 0.275 | 0.003 |
+------+-------+-------+
####
use PDL::Apply 'apply_rolling';
sub sms_WxH_PDL_ar ( $m, $w, $h ) {
$m-> apply_rolling( $w, 'sum' )
-> slice( [ $w - 1, -1 ], [] )
-> xchg( 0, 1 )
-> apply_rolling( $h, 'sum' )
-> slice( [ $h - 1, -1 ], [] )
-> xchg( 0, 1 )
}
my $m = random 150, 150;
$t = time;
sms_WxH_PDL_ar( $m, 2, 2 );
say time - $t;
$t = time;
sms_WxH_PDL_ar( $m, 10, 10 );
say time - $t;
$t = time;
sms_WxH_PDL_ar( $m, 75, 75 );
say time - $t;
__END__
3.36784410476685
3.20836806297302
1.28170394897461