Automation of for loops and array dimesion

5 views (last 30 days)
Hello members, I need some help.
rho_grid=[-1,0,1];
Nrho=length(rho_grid);
for i=Nrho:-1:1
rho1=rho_grid(i);
for j=Nrho:-1:1
rho2=rho_grid(j);
[U,S,V]= svd(pro(rho1,rho2));
[Uvals(:,:,j,i),Svals(:,:,j,i),Vvals(:,:,j,i)]= svd(pro(rho1,rho2)); % I need this later for plotting
Ti= inv(sqrt(S))*U'*Rq(rho1,rho2);
T = Rp(rho1,rho2)*V*inv(sqrt(S));
systrans(:,:,j,i)= ss2ss(sys.Data(:,:,j,i),Ti);
end
end
1) I have 2 'for' loops now to determine value of parameters 'rho1' and 'rho2'. In my case, this parameters will increase to very large number. If I have rho1,....,rho10 , I dont want to create all extra for loops manually. Is there any way that I can automate it? Please note that rho_grid = [-1,0,1] is same for all parameters. Similarly,
2) The dimension of 'systrans' will also change. If I add 1 extra for loop to this I have to change from my current systrans(:,:,j,i) to
systrans(:,:,k,j,i)
which I want to avoid and automate.
3) Last question, I want plot values stored in 'Svals' over parameter range. When I had only one parameter 'rho1' I did this:
plot(rho_grid,reshape(Svals(i,i,:),[1,length(rho_grid)]));
which gave me 2D plot. for 2 parameters case, I want to get a 3D plot and see the variation of entries of Svals over rho1 and rho2. How can I do that? I have attached Uvals,Svals and Vvals in .mat file
  5 Comments
Sandeep Parameshwara
Sandeep Parameshwara on 23 Jan 2020
Hi Mr. Guillaume, you are right. pro is just a double matrix of size 8*8 in above case. Even though I increase number of parameters, it will be 2D matrix. I wrote Uvals(:,:,j,i) etc to store the values for each value of rho as you can see from my code.

Sign in to comment.

Accepted Answer

Guillaume
Guillaume on 24 Jan 2020
Assumption: the matrices returned by pro, Rq, Rp, and ss2ss are all the same size (size defined by pro_outsize below
I've not changed your equations. As mentioned before inv is not recommended.
Untested code. I may have made a mistake but the princiiple should be correct:
num_rhos = 10; %number of rho variables
rho_grid = [-1 0 1];
pro_outsize = [8 8]; %size of matrices returned by pro, Rq, Rp and ss2ss
%create all permutation of picking num_rhos elements out of rho_grid. There's numel(rho_grid) ^ num_rhos permutations
rhoperms = cell(1, num_rhos);
[rhoperms{:}] = ndgrid(rho_grid);
rhoperms = reshape(cat(num_rhos + 1, rhoperms{:}), [], num_rhos); %each row of rows is a unique permutation. Height of the matrix is numel(rho_grid) ^ num_rhos
%preallocation of results as 3D matrices. Can be reshape to ND matrices at the end
U_vals = zeros([pro_outsize, size(rhoperms, 1)]);
S_vals = U_vals;
V_vals = U_vals;
systrans = U_vals;
%reshape sys.Data into 3D matrix for easier indexing later on
sysdata = reshape(sys.Data, size(sys.Data, 1), size(sys.Data, 2), []);
%loop over each unique combination
for ridx = 1:size(rhoperms, 1)
rhos = num2cell(rhoperms(ridx, :));
[U_vals(:, :, ridx), S_vals(:, :, ridx), V_vals(:, :, ridx)] = svd(pro(rhos{:}));
Ti = inv(sqrt(S_vals(:, :, ridx))) * U_vals(:, :, ridx)' * Rq(rhos{:}); %probably, shouldn't use inv
T = Rp(rhos{:}) * V_vals(:, :, ridx) * inv(sqrt(S_vals(:, :, ridx)));
systrans(:, :, ridx) = ss2ss(sysdata(:, :, ridx), Ti);
end
%Optionally: reshape into N-D matrices, where N = num_rhos + 2
newshape = [pro_outsize, repelem(numel(rho_grid), num_rhos)];
U_vals = reshape(U_vals, newshape);
S_vals = reshape(S_vals, newshape);
V_vals = reshape(V_vals, newshape);
systrans = reshape(systrans, newshape);
The code use the expansion of cell arrays into comma-separated lists to build the inputs of pro, Rq, and Rp.
  3 Comments
Guillaume
Guillaume on 24 Jan 2020
I'm not sure I completely follow, since the rhos are the same in any dimension, wouldn't (:, :, 1, 3) and (:, :, 3, 1) be the same?
You can reorder the dimensions with permute, so after the reshape, you can do:
systrans = permute(systrans, [1 2, num_rhos+2:-1:2]);
I think! It's always difficult to think in >3D. But again, I'm not sure it makes a difference.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by