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 ```