Looping through a matrix of unknown dimensions
52 views (last 30 days)
Show older comments
Hi everyone, I know similar questions have been asked before but I can't find anything that would be relevant for my particular situation. I need to write a function that would retrieve a vector from a matrix of unknown size, operate on it, and then return the results into an output matrix of size equal to the starting matrix. Here is an example of what I need:
start_matrix = rand(5,6,7,8,1000);
s = size(start_matrix);
output_matrix = zeros(s);
for a = 1:s(1)
for b = 1:s(2)
for c = 1:s(3)
for d = 1:s(4)
operate_data = start_matrix(a,b,c,d,:); % retrieving the vector of length=1000
outputmatrix(a,b,c,d,:) = fft(operate_data); % using fast fourier transform as an example, there are a few more calculations here
end
end
end
end
Now the number of dimensions can vary and the size of them can vary too, meaning that the number of loops needed will be variable. So my question is, how do I make this work with e.g. size = [10,15,20,500], where I loop through the first three dimensions, doing calculations on vectors of length = 500?
I've been trying to implement something like in this example by Jan:
But I don't know how to make it work when I am not operating on each element, but rather a vector, each time.
Help is greatly appreciated! Thanks.
Accepted Answer
Jan
on 28 Feb 2019
Edited: Jan
on 28 Feb 2019
What about reshaping the array?
start_matrix = rand(5,6,7,8,1000);
s = size(start_matrix);
M = reshape(start_matrix, [], s(end));
Mfft = fft(M, [], 2);
Result = reshape(Mfft, s);
isequal(Result, output_matrix)
>> logical(1)
This is much faster than a loop. But if you really need a set of loop with distinct dimensions, this is the way to go:
X = rand(5,6,7,8,1000);
Out = zeros(size(X)); % Pre-allocate the output
% Prepare index vector as cell, e.g. {1,2,3,4,':'} :
nv = ndims(X) - 1; % Exclude last dimension
v = [repmat({1}, 1, nv), {':'}];
vLim = size(X); % Loop until size (last element not used)
ready = false;
while ~ready
% The actual operation:
Out(v{:}) = fft(X(v{:}));
% Update the index vector:
ready = true; % Assume that the WHILE loop is ready
for k = 1:nv % Loop through dimensions
v{k} = v{k} + 1; % Increase current index by 1
if v{k} <= vLim(k) % Does it exceed the length of current dim?
ready = false; % No, WHILE loop is not ready now
break; % v(k) increased successfully, leave "for k" loop
end
v{k} = 1; % Reset v{k}, proceed to next k
end
end
The idea is to create a cell vector, which contains the indices, e.g. {1,2,3,4,':'}. Then this is equivalent:
index = {1, 2, 3, 4, ':'}
A(1,2,3,4,:)
A(index{:})
Now tis index vector is used and updated in the for loop: The first element is increased by 1. If it does not exceed size(X, 1), the next iteration is performed. If it exceeds size(X, 1), the element is reset to 1 and the next one is increased.
3 Comments
More Answers (2)
Guillaume
on 28 Feb 2019
Edited: Guillaume
on 28 Feb 2019
I can see two ways of solving your problem.
The first way, would be to permute the dimensions, so that your : is move to the first dimension, This way, you can linear indices. For m a matrix of size MxNxPx..., m(:, 1, 1, 1, ...) is the same as m(1:M), m(:, 2, 1, 1, ...) is m((1:M) + M), etc. so:
m = rand(5,6,7,8,1000);
%move the colon dimension to the first one:
m = permute(m, [ndims(m), 1:ndims(m) - 1]);
output = zeros(size(m));
%use linear indexing
sz = size(m);
indices = 1:sz(1);
for stride = 1:prod(sz(2:end))
operate_indices = indices + (stride-1)*sz(1);
output(operate_indices) = fft(m(operate_indices));
end
%at the end you can restore the original ordering if you want:
output = permute(output, [2:ndims(m), 1]);
The other way is to use the expansion of cell arrays into comma-separated list and store the indices into a cell array:
m = rand(5,6,7,8,1000);
%calculate all indices excluding last dimension
sz = size(m)
indices = arrayfun(@(s) 1:s, sz(1:end-1), 'UniformOutput', false);
[indices{:}] = ndgrid(indices{:}); %cartesian product of all indices of all dimensions. Using expansion of cell array into csl
%convert into a 2D matrix
indices = reshape(cat(numel(indices)+1, indices{:}), numel(indices{1}), [];
% at this point each row of indices, is a list of index of each dimension. Split back into cell array
indices = num2cell(indices);
%iterate over the indices combination
output = zeros(size(m));
for row = 1:size(indices, 1)
output(indices{row, :}, :) = fft(m(indices{row, :}, :));
end
Andrei Bobrov
on 28 Feb 2019
For your case:
start_matrix = rand(5,6,7,8,1000);
nd = 1:ndims(start_matrix);
a = permute(start_matrix,nd([end,1:end-1]));
out_matrix = permute(fft(a),nd([2:end,1]));
0 Comments
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!