Automation of for loops and array dimesion

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

[Guillaume]

It's really not my domain, so I can't really tell what you're doing but note that computing inverse is rarely the way to go, see the tips section in the documention of inv.

With regards to your question, does this mean that the pro function (I assume it's a function) will then take 10 inputs (not a good idea) or a vector of 10 elements (much better)? Will the subsequent code be designed to work with 10-D systrans, Uvals, Svals, and Vvals? Can the function ss2ss cope with a 10-D matrix?

[Sandeep Parameshwara]

Hi,thanks for the reply. I will look into the inverse and change it to backslash. Regarding your other queries, 1)pro is a function handle,so i need to provide argument list. So I have to give a list of rho1,..,rho10 . For time being I am using 'eval' (pls don't judge me,I know that it's not recommended. I am planning to switch to more sophisticated method soon) 2) yes it will be 10D systrans,Uvals,Svals and Vvals, ss2ss can cope that. I am not planning to go higher than 10, but i want to automate this code nevertheless. I need to switch between number of parameters often,so it's impractical for me to edit it everytime. Greatful to get any help, hints from you.

[Mohammad Sami]

Perhaps you can do it recursively. Define a base case of 4-D array and if its higher, it can call itself with 1 less dimension.

[Guillaume]

I wouldn't use a recursion, it can all be done with a single loop. One last question, what's the size of the matrix returned by pro (or what's the size of U, S, V)? I'm assuming that it's independent of the number of rhos but if it does depend what's the formula.

It's required in order to preallocate, U_vals, S_sals, V_vals and systrans

[Sandeep Parameshwara]

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.

Answer 1

[Guillaume]

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.

[Sandeep Parameshwara]

Hi,

Thanks a lot for this help. It almost works. one last question, there is one problem. sys.Data is 4D data. Its '4D ss' (ss= state space) data. Hence when I reshape it it becomes 5D ss data. Hence, I decided not to do it and used sys.Data in ss2ss command. However after that I use a command in my code which takes systrans data and maps it to grid. The actual command is

sysbal=pss(systrans,rho_rgrid)

where rho_rgrid is given by the command

rho_rgrid=rgrid(rho{1},rho{2});

for 2 parameters. These are not matlab command but from external toolbox. I think below figure explains the problem. systrans is in the form systrans(:,:,1,1),systrans(:,:,2,1),...,systrans(:,:,9,1) where as I need systrans(:,:,1,1),systrans(:,:,1,2),systrans(:,:,1,3),systrans(:,:,2,1),...,systrans(:,:,3,3)

so sysbal=pss(systrans,rho_rgrid) fails and gives me the error

%Error using pss (line 71)
%Dimensions of Domain and Data are incompatible.

domain menas rho_rgrid and data is systrans. So could you please tell me how I can map it to the grid?

Thanks a lot

[Sandeep Parameshwara]

I reshaped it using reshape(systrans,[3 3]), it works. How should I make sure that the reshaped sequence is as intended? Sorry for repeated questions.

[Guillaume]

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.