Strategies to optimize nested for loops?

I have the following code in which I have to loop through multiple indices. Each for loop has a length of 3, 363, 10, 100, and 11. I'm doing this because I'm performing a Monte Carlo simulation on multiple functions of random variables and want to retain the looping calculation. This obviously takes quite a bit of time. Any suggestions for strategies around nested for loops?

% MCS for CPT of Flow Volume Node and Cumulative Nodes
for iDM=1:length(DamState)  % Number of flow area states
    for iCombFGH=1:size(IndMatFGH,1) % External flood height states
        for iCombV1 = 1:length(VBinEdges_Md) % Single flow volume state
            for iCombCum=1:size(IndMatVCum,1) % Cumulative flow volume state combos
                for iCombIFH=1:length(IBinEdges_Md) % States of IFH
                
                % Flow area sim, 3 damage states
                Sig=[]; Sig=0.2*Amean;
                pd_Asim=[]; pd_Asim = makedist('normal',Amean(iDM),Sig(iDM));
                Asim = []; Asim=random(pd_Asim,[Nsim,1]);
                Asim(Asim<0)=0; % make sure the area is postive
                Asim(Asim>Aopen)=Aopen;  %make sure area is never bigger than the full area
    
                % External Flood Height sim, 11 states
                a=[];
                c=[];
                pdEsim=[];
                Esim=[];
                % Define distribution of external flood height within bin
                a=ExtSurgeBinEdges_Lo(IndMatFGH(iCombFGH,3));
                b=ExtSurgeBinEdges_Md(IndMatFGH(iCombFGH,3));
                c=ExtSurgeBinEdges_Hi(IndMatFGH(iCombFGH,3));
                pdEsim=makedist('Uniform',a,c);
                Esim=random(pdEsim,[Nsim,1]); %simulate external flood heights within bin (Uniform)
                Esim=Esim.*(Esim>=n_ins); % make sure the difference is never below installation heigh
    
                % Flowrate calculation
                Qsim=[]; Qsim=Cd.*Asim.*sqrt(2*g*(Esim)); %ft^3/s
                
                % Flood duration sim
                % [Note: Currently distribution is fixed; eventually should be a function of flood height] 
                Mu=[]; Mu=3600; %sec
                Sig=[]; Sig=0.2*Mu;
                pdDsim=[]; pdDsim=makedist('normal',Mu,Sig); 
                Dsim=[]; Dsim=random(pdDsim,[Nsim,1]); %simulate duration
                
                % Volumetric inflow calculation
                Vsim=[]; Vsim=Qsim.*Dsim; %ft^3
                Vsim(Vsim<0)=0; %Ensure positive values
                Vsim=Vsim.*(Vsim<=Vroom(:,1)); %Ensure Vsim is less than Vroom (3 Vroom configs)
                
                % PMF_Vol for each floor elevation configuration
                %Floor at plant elevation
                Hist=[]; [Hist,~] = histcounts(Vsim(:,1),VBinEdges);
                Hist=Hist';
                PMF_V_1(:,iDM)=Hist./sum(Hist);
                
                % Cum Flow 1 calculation
                Vcum1sim = Vsim(:,1)+Vsim(:,1);
                Vcum1sim(Vcum1sim>Vroom(:,1))=Vroom(:,1);  %make sure volume is never bigger than the room volume
            
                Hist=[]; [Hist,~] = histcounts(Vcum1sim,VBinEdges);
                Hist=Hist';
                PMF_cumV1(:,iCombCum)=Hist./sum(Hist);
                
                % For cum flow 2 calculation (three sources)
                Vcum2sim = Vcum1sim+Vsim(:,1);
                Vcum2sim(Vcum2sim<0)=0; % make sure the area is postive
                Vcum2sim(Vcum2sim>Vroom(:,1))=Vroom(:,1);  %make sure volume is never bigger than the room volume
            
                Hist=[]; [Hist,~] = histcounts(Vcum2sim,VBinEdges);
                Hist=Hist';
                PMF_cumV2(:,iCombCum)=Hist./sum(Hist);
                
                % For cumflow 3 (four sources)
                Vcum3sim = Vcum2sim+Vsim(:,1);
                Vcum3sim(Vcum3sim<0)=0; % make sure the area is postive
                Vcum3sim(Vcum3sim>Vroom(:,1))=Vroom(:,1);  %make sure area is never bigger than the room volume
            
                Hist=[]; [Hist,~] = histcounts(Vcum3sim,VBinEdges);
                Hist=Hist';
                PMF_cumV3(:,iCombCum)=Hist./sum(Hist);

                % Internal flood height calculation
                IFHsim = Vcum3sim./(int_W*int_L);

                Hist=[]; [Hist,~] = histcounts(IFHsim,IBinEdges);
                Hist=Hist';
                PMF_cumV3(:,iCombIFH)=Hist./sum(Hist);
                end 
            end 
        end
    end
end 

Answer 1

[Daniel Bengtson]

If you are going to be looping many times through a set of operations with a goal of execution time or efficiency, you should make sure that only those operations that are absolutely neccessary are in the loop.

It appears that there are sections of that code that do not change loop by loop. Is there any reason that this part even has to be in the loop at all rather than just being run once at the start?

                % Flood duration sim
                % [Note: Currently distribution is fixed; eventually should be a function of flood height] 
                Mu=[]; Mu=3600; %sec
                Sig=[]; Sig=0.2*Mu;
                pdDsim=[]; pdDsim=makedist('normal',Mu,Sig); 
                Dsim=[]; Dsim=random(pdDsim,[Nsim,1]); %simulate duration

Can these three functions be combined into one function that returns 3 outputs? I assume at least some of the operations being done in them overlap, so it may allow for some time savings to combine them and make use of the operations you have already done to speed up later ones.

                a=ExtSurgeBinEdges_Lo(IndMatFGH(iCombFGH,3));
                b=ExtSurgeBinEdges_Md(IndMatFGH(iCombFGH,3));
                c=ExtSurgeBinEdges_Hi(IndMatFGH(iCombFGH,3));

This whole section looks like it could be moved up to only run on change of 'iCombFGH'

                                % Flow area sim, 3 damage states
                Sig=[]; Sig=0.2*Amean;
                pd_Asim=[]; pd_Asim = makedist('normal',Amean(iDM),Sig(iDM));
                Asim = []; Asim=random(pd_Asim,[Nsim,1]);
                Asim(Asim<0)=0; % make sure the area is postive
                Asim(Asim>Aopen)=Aopen;  %make sure area is never bigger than the full area
    
                % External Flood Height sim, 11 states
                a=[];
                c=[];
                pdEsim=[];
                Esim=[];
                % Define distribution of external flood height within bin
                a=ExtSurgeBinEdges_Lo(IndMatFGH(iCombFGH,3));
                b=ExtSurgeBinEdges_Md(IndMatFGH(iCombFGH,3));
                c=ExtSurgeBinEdges_Hi(IndMatFGH(iCombFGH,3));
                pdEsim=makedist('Uniform',a,c);
                Esim=random(pdEsim,[Nsim,1]); %simulate external flood heights within bin (Uniform)
                Esim=Esim.*(Esim>=n_ins); % make sure the difference is never below installation heigh
    
                % Flowrate calculation
                Qsim=[]; Qsim=Cd.*Asim.*sqrt(2*g*(Esim)); %ft^3/s
                
                % Flood duration sim
                % [Note: Currently distribution is fixed; eventually should be a function of flood height] 
                Mu=[]; Mu=3600; %sec
                Sig=[]; Sig=0.2*Mu;
                pdDsim=[]; pdDsim=makedist('normal',Mu,Sig); 
                Dsim=[]; Dsim=random(pdDsim,[Nsim,1]); %simulate duration
                
                % Volumetric inflow calculation
                Vsim=[]; Vsim=Qsim.*Dsim; %ft^3
                Vsim(Vsim<0)=0; %Ensure positive values
                Vsim=Vsim.*(Vsim<=Vroom(:,1)); %Ensure Vsim is less than Vroom (3 Vroom configs)

These are just what popped out at first glance. A good first step in optimizing would be to run the code profiler and use the flame graph to identify where the most time is being spent. Attacking the areas that either consume the most time, or are run the most number of times will likely lead to the easiest performance increases.

[Jeff Miller]

Similarly, these two lines:

                Sig=[]; Sig=0.2*Amean;
                pd_Asim=[]; pd_Asim = makedist('normal',Amean(iDM),Sig(iDM));

can follow immediately after the 'for iDM' loop. They produce the same result across all combinations of the other indices, so there is no need to execute them repeatedly. You will need the new random numbers for each combination, though.

You might also look at parfor--this looks like a good candidate for that if you have multiple processors.

[Joy]

@Daniel Bengtson @Jeff Miller

Thanks for the suggestions, I'm looking into the code profiler and parfor to see if those help.

Jeff, what do you mean by I'll need new random numbers for each combination when I move the code right underneath their associated for loops?

Thanks

[Jeff Miller]

@Joy

I just meant that the statement

Asim=random(pd_Asim,[Nsim,1]);

must still stay in the inner-most for loop even though you move the pd_Asim makedist command to an outer level of the nested for's.