You are looping because you are not trusting the threading engine.
The strategy is to push your unused dimensions above the active
dimensions of your operator, and then let them ride along. Matrix
multiplication has two active dimensions, so you have to add dummy
dimensions to make your data into a collection of (1,n,x,y) column
vectors and (n,1,x,y) row vectors. Then you can just carry out your
multiplication. You probably want to get out of the habit of using
'slice' and into the habit of using "PDL::NiceSlice", which yields a
much cleaner syntax.
$colvecs = $allims->(*1); # colvecs dims: (1,n,x,y)
$rowvecs = $allims->(:,*1); # rowvecs dims: (n,1,x,y)
$result = $rowvecs x $matrix x $colvecs; # result dims: (1,1,x,y)
$out = $result->clump(3); # out dims: (x,y)
A little more concise:
$vecs = $allims->(*1);
$out = ($vecs->transpose x $matrix x $vecs)->clump(3);
You can even do it in one expression if you prefer:
$out = ($allims->(:,*1) x $matrix x $allims->(*1))->clump(3);
BTW, your '->cat' call is doing nothing, since it sticks together a
bunch of PDLs and you're just giving it one PDL.
Cheers,
Craig
On Apr 15, 2008, at 12:49 PM, Ashish Mahabal wrote:
>
> Hi,
>
> I know there has to be a better, more succinct way of doing this,
> but I
> just don't seem to get it, not having done any explicit threading:
>
> $allims is a stack of n images of size (xsize,ysize) from which I
> need to
> pull out vectors that have the same pixel for each image in the
> stack and
> do some matrix operations on it. The loops clearly slow it down. How
> do I
> get rid of the loops? ($invcovari is a nxn matrix).
>
> my $out= ones($xsize, $ysize);
> for my $j (0..$xsize-1){
> for my $k (0..$ysize-1){
> my $vecf = (slice $allims,"*,$j,$k")->flat->cat;
> my $vecft = $vecf->transpose;
> (slice $out,"$j,$k") .= $vecf x ($invcovari x $vecft);
> }
> }
>
> Thanks.
>
> -ashish
>
> Ashish Mahabal, Caltech Astronomy, Pasadena, CA 91125
> http://www.astro.caltech.edu/~aam aam at astro.caltech.edu
>
> Why is "abbreviation" such a long word?
>
> _______________________________________________
> Perldl mailing list
> [email protected]
> http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
>
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
