Thread Subject: multidimensional indexing

Hi,

I'm looking for a sort of generalization of diag for
multidimensional arrays. Specifically suppose I have:

A(:,:,1) = [1 2 3; 1 2 3; 1 2 3];

A(:,:,2) = [4 5 6; 4 5 6; 4 5 6];

A(:,:,3) = [7 8 9; 7 8 9; 7 8 9];

Is there a quick way to get

B = [1 5 9; 1 5 9; 1 5 9], so the first column of B is the
first column of the first page, the second column is the
second column of the second page, etc? (avoiding for
loops) Thanks!

"Darius " <martin@econ.ucsb.edu> wrote in message
<g0b6eb$18j$1@fred.mathworks.com>...
> Hi,
>
> I'm looking for a sort of generalization of diag for
> multidimensional arrays. Specifically suppose I have:
>
> A(:,:,1) = [1 2 3; 1 2 3; 1 2 3];
>
> A(:,:,2) = [4 5 6; 4 5 6; 4 5 6];
>
> A(:,:,3) = [7 8 9; 7 8 9; 7 8 9];
>
> Is there a quick way to get
>
> B = [1 5 9; 1 5 9; 1 5 9], so the first column of B is
the
> first column of the first page, the second column is the
> second column of the second page, etc? (avoiding for
> loops) Thanks!
>

One of many methods:

i = repmat((1:size(A,1)),size(A,1),1);
I = sub2ind(size(A),i.',i,i);
A(I)

"Darius ":
<SNIP repeating his/her ND-diag...

one of the many solutions

% the data
     m=reshape(1:27,[3,3,3])
% the engine
     r=ones(3,1)*arrayfun(@(x) m(x,x,x),1:3)
%{
     r=
     1 14 27
     1 14 27
     1 14 27
%}

us

"us " <us@neurol.unizh.ch> wrote in message <g0bigi$cp3
$1@fred.mathworks.com>...
> "Darius ":
> <SNIP repeating his/her ND-diag...
>
> one of the many solutions
>
> % the data
> m=reshape(1:27,[3,3,3])
> % the engine
> r=ones(3,1)*arrayfun(@(x) m(x,x,x),1:3)
> %{
> r=
> 1 14 27
> 1 14 27
> 1 14 27
> %}
>
> us

Ah but this doesn't quite satisfy his request which is the
equivalent to:

r = [m(:,1,1) m(:,2,2) m(:,3,3)];

This would require the following change to your code:

r=cell2mat(arrayfun(@(x) m(:,x,x),1:3,'uni',false));

And, of course the comparison between the methods (and an
additional one using FOR-loops):

% The data
N = 100;
m = rand(N,N,N);
t = {'helper'; 'us'; 'FOR-loops'};

% The tests
tic
for n = 1:1000
i = repmat((1:size(m,1)),size(m,1),1);
I = sub2ind(size(m),i.',i,i);
r = m(I);
end
t{1,2} = toc;

tic
for n = 1:1000
r = cell2mat(arrayfun(@(x) m(:,x,x),1:3,'uni',false));
end
t{2,2} = toc;

tic
for n = 1:1000
  r = zeros(N,N);
  for j = 1:N
    r(:,j) = m(:,j,j);
  end
end
t{3,2} = toc;

disp(t)

    'helper' [1.1061]
    'us' [0.5106]
    'FOR-loops' [0.2666] % The winner


FOR-loops win again!

"helper ":
<SNIP scolding <us>'s solution... :-)

> Ah but this doesn't quite satisfy his request which is
the equivalent to:
> r = [m(:,1,1) m(:,2,2) m(:,3,3)];

<helper>, i couldn't agree more! - but then: i looked at
the OP's example, in particular his/her (ambiguous?)
solution

     B = [1 5 9; 1 5 9; 1 5 9]
%{
     1 5 9
     1 5 9
     1 5 9
%}

that had me come up with the <repeating his/her ND-
diag...>...
but again, you're most likely correct! - in which case
another solution might look like this

% the data
     m=repmat(magic(3),[1,1,3]);
% the engine
     r=reshape(m,size(m,1),[]);
     r=r(:,1:size(m,2)+1:end)
%{
     8 1 6
     3 5 7
     4 9 2
%}

now, using your nice timing snippet (below), things look
like this on a wintel: ic2/2*2.4mhz/2g/winxp.sp2/r2007b)

     'helper' [0.77992]
     'us' [0.054153]
     'FOR-loops' [0.21931]

again, just for fun :-)
us

% The data
     N = 100;
     m = rand(N,N,N);
     t = {'helper'; 'us'; 'FOR-loops'};

% helper
     tic
for n = 1:1000
     i = repmat((1:size(m,1)),size(m,1),1);
     I = sub2ind(size(m),i.',i,i);
     r = m(I);
end
     t{1,2} = toc;
     r1=r;
% us
     tic
for n = 1:1000
     r=reshape(m,size(m,1),[]);
     r=r(:,1:size(m,2)+1:end);
end
     t{2,2} = toc;
     r2=r;
% loop
     tic
for n = 1:1000
     r = zeros(N,N);
for j = 1:N
     r(:,j) = m(:,j,j);
end
end
     t{3,2} = toc;
     r3=r;
isequal(r1,r2,r3)
% ans = 1
     disp(t)

Thank you all! There's still so much to learn.

"us " <us@neurol.unizh.ch> wrote in message
<g0c3nr$78u$1@fred.mathworks.com>...
> "helper ":
> <SNIP scolding <us>'s solution... :-)
>
> > Ah but this doesn't quite satisfy his request which is
> the equivalent to:
> > r = [m(:,1,1) m(:,2,2) m(:,3,3)];
>
> <helper>, i couldn't agree more! - but then: i looked at
> the OP's example, in particular his/her (ambiguous?)
> solution
>
> B = [1 5 9; 1 5 9; 1 5 9]
> %{
> 1 5 9
> 1 5 9
> 1 5 9
> %}
>
> that had me come up with the <repeating his/her ND-
> diag...>...
> but again, you're most likely correct! - in which case
> another solution might look like this
>
> % the data
> m=repmat(magic(3),[1,1,3]);
> % the engine
> r=reshape(m,size(m,1),[]);
> r=r(:,1:size(m,2)+1:end)
> %{
> 8 1 6
> 3 5 7
> 4 9 2
> %}
>
> now, using your nice timing snippet (below), things look
> like this on a wintel: ic2/2*2.4mhz/2g/winxp.sp2/r2007b)
>
> 'helper' [0.77992]
> 'us' [0.054153]
> 'FOR-loops' [0.21931]
>
> again, just for fun :-)
> us
>
> % The data
> N = 100;
> m = rand(N,N,N);
> t = {'helper'; 'us'; 'FOR-loops'};
>
> % helper
> tic
> for n = 1:1000
> i = repmat((1:size(m,1)),size(m,1),1);
> I = sub2ind(size(m),i.',i,i);
> r = m(I);
> end
> t{1,2} = toc;
> r1=r;
> % us
> tic
> for n = 1:1000
> r=reshape(m,size(m,1),[]);
> r=r(:,1:size(m,2)+1:end);
> end
> t{2,2} = toc;
> r2=r;
> % loop
> tic
> for n = 1:1000
> r = zeros(N,N);
> for j = 1:N
> r(:,j) = m(:,j,j);
> end
> end
> t{3,2} = toc;
> r3=r;
> isequal(r1,r2,r3)
> % ans = 1
> disp(t)