Path: news.mathworks.com!newsfeed-00.mathworks.com!panix!bloom-beacon.mit.edu!llnews!53ab2750!not-for-mail
From: Peter Boettcher <boettcher@ll.mit.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: vectorise a reshape loop
References: <g5krk2$77v$1@fred.mathworks.com>
Message-ID: <muy8wvv8dy0.fsf@G99-Boettcher.llan.ll.mit.edu>
Organization: MIT Lincoln Laboratory
User-Agent: Gnus/5.110006 (No Gnus v0.6) Emacs/23.0.0 (gnu/linux)
Cancel-Lock: sha1:o6uNW1JD0DjUHZzPF0NHnb4v1vY=
MIME-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Lines: 115
Date: Mon, 21 Jul 2008 11:16:55 -0400
NNTP-Posting-Host: 155.34.163.114
X-Complaints-To: news@ll.mit.edu
X-Trace: llnews 1216652673 155.34.163.114 (Mon, 21 Jul 2008 11:04:33 EDT)
NNTP-Posting-Date: Mon, 21 Jul 2008 11:04:33 EDT
Xref: news.mathworks.com comp.soft-sys.matlab:480685



"Dave Brackett" <davebrackett@hotmail.com> writes:

>> The size of the left side IS NOT 2x2x90x1. q is a scalar,
>> so the left side is 2x2x1x1. The right side, I assume, is
>> 2x2x1x16. But don't guess or believe me. If you can't 
>> figure out your dimension issues, I meant what I said 
>> about "if you're not sure, assign it to a variable and 
>> ask MATLAB to tell you what the size is".
>> tmp = ndfun('mprod', M(:,:,q,:));
>> size(tmp)
>
> I am getting quite confused with what I thought would be a
> simple thing. prodH and M are contained within a loop of 90
> iterations, so the size of the resulting prodH matrix after
> running the loop should be 2x2x90x1.
>
> M is of the size 2x2x90x16 and i am iterating through this
> using a loop to find the product of the 4th dimension
> (length 16) for each of the indices in the 3rd dimension
> (length 90). 
>
> This would then result in prodH of the size 2x2x90x1 with
> each 2x2 matrix being the product of the 16 in M.
>
> I hope that makes sense. That is what I am trying to achieve
> and thought that ndfun('mprod') would be enable that. If my
> loop is confusing and seems to suggest something different
> then apologies.
>
> Please could you show me how to write this correctly as I am
> starting to go mad with it? Thanks.

Until now, I have been unable to understand exactly what you want where.
With this explanation, I think I can help.

I admit that I am still frustrated that you refuse to follow my
debugging advice by looking at the size of the matrix generated by
ndfun, and then reading the help for ndfun.  So instead of just writing
the answer, I'll go through the exercise here, hopefully providing an
illustration of the debugging process I was asking you to try.

M = rand(2,2,90,16);
tmp = ndfun('mprod', M(:,:,1,:));
size(tmp)

ans = 
   2  2  16


Hmm, why is it 2x2x16?

"help ndfun"
    'mprod' behaves differently.  It cumulatively multiplies a set
    of matrices and produces a single 2D output.  The equivalent
    code is:
        C = A(:,:,1);
        for i=2:N
            C = C * A(:,:,i);
        end
    2D inputs return themselves.  Inputs with more than 3
    dimensions collapse the 3rd dimension only.  So with an A of
    size [2 2 7 3 4],
        C = ndfun('mprod', A);
    is equivalent to
        for i=1:3
            for j=1:4
                C(:,:,i,j)=ndfun('mprod',M(:,:,:,i,j));
            end
        end
    and C will have size [2 2 3 4].


Reading that carefully, especially the second example, suggests that
inputting a 2x2x1x16 matrix will produce a 2x2x16 output, essentially
doing nothing.  Since ndfun('mprod') works down the third dimension, you
want your x16 size to be on that third dimension.  "Squeeze" will do
that in this case, but it's just a special case of permute, so whichever
you want.

So this should work:

prodH=zeros(2,2,length(freq),1);
for p=1:length(freq)
    q=p;
    prodH(:,:,q,1)=ndfun('mprod', squeeze(M(:,:,q,:)))
end

Let's check:

Right side:

M(:,:,q,:) is 2x2x1x16
squeeze(M(:,:,q,:)) is 2x2x16
   (you can always ignore trailing 1 dimensions)
output of ndfun is 2x2

Left side:
prodH(:,:,q,1) is 2x2

Good!

OK, one more thing.  Again looking carefully at the help for ndfun, the
example shows a case where the mprod is repeated for higher-dimension
indices.  So it appears that if you put your x90 dimension at the end
(in the fourth dimension), ndfun will automatically do what you want,
which is to collapse the 3rd dimension (16) by multiplying the 2x2
matrices, then do this for all of the 4d "pages".  So you want to switch
3rd and 4th dimensions of M:

This should do everything, sans loop:

prodH = ndfun('mprod', permute(M, [1 2 4 3]));


-Peter