I have the following code which is too slow right now and I am looking for better ways to write it.
This code is part of a much longer code including a for loop t=1:T, therefore the following code runs for every t.
I have F firms, which produce Yd. In every period t they buy a certain amount K of machine tools, whose productivity is given by A; therefore, the productive capacity of the machine tools bought in t equals AK=A.*K. Firms' total productive capacity equals sum(AK). Since sum(AK) might be greater than what they want to produce Yd, firms have to determine which fraction they will utilize (omega) of each machine "v" they own, considering that they start using the most productive machine tools first (those with the highest A). Thus, starting from the most productive machines, if the productive capacity from a certain machine is lower than Yd, omega of that machine will be set to 1, if the productive capacity has already reached Yd, omega of that machine will be set to 0, otherwise omega will be set to the fraction required. How can I code that?
This is what I did but it is definitely too slow:
T=100;
F=50;
A=rand(T,F); %initialized randomly (not as in my real code)
K=rand(T,F); %idem
AK=A.*K; %productivity of machine tools
Yd=rand(1,F); %initialized randomly
omega=zeros(T,F); %needs to be computed
OAK(:,:)=zeros(T,F); %needs to be computed
for f=1:F
[A_sorted_value,A_sorted_index]=sort(A(:,f)); %sort machines according to their productivity level
while sum(A_sorted_value)>0 && sum(OAK(1:t,f))<Yd(f) %apply rule as long as there are machines available and desired production has not been achieved
v=A_sorted_index(end); %pick the best machine of capital, i.e. highest productivity A
if AK(v,f)+sum(OAK(1:t,f))<=Yd(f) %assign omega=1 if existing capacity is below desired production
omega(v,f)=1;
OAK(v,f)=omega(v,f).*AK(v,f); %update existing capacity
elseif AK(v,f)+sum(OAK(1:t,f))>Yd(f) && Yd(f)>sum(OAK(1:t,f)) %assign omega so that machine v can fill desired production
omega(v,f)=(Yd(f)-sum(OAK(1:t,f)))/AK(v,f);
OAK(v,f)=omega(v,f).*AK(v,f);
end
A_sorted_index(end)=[]; %remove machine v from the list and pick the next one if condition "while" is satisfied
A_sorted_value(end)=[]; %otherwise go to the next firm
end
end
I hope it is clear and thanks in advance!
It appears that you are running some kind of simulation. This means (if results depend on previous iterations), your t and f loop can't be optimized away. So the only optimization target available would be the inner while loop. I'd recommend you consider whether or not you can pre-calculate or vectorize anything (check your math on paper).
But since I can't really help you with that, here are some things which might help.
One low hanging fruit I can offer immediately, is you can use the sort function to sort the entire matrix in some set direction, so you can move that outside the f loop.
% Sort the entire A array, column-wise.
[AVals, AIdxs] = sort(A,1);
for f = 1:F
A_sorted_value = AVal(:,f);
A_sorted_index = AIdx(:,f);
...
end
Next, I recommend that you use a for loop instead of your while loop.
Since you already know that you will be iterating through the entire sorted A list, or until another condition (sum(OAK(1:t,f))<Yd(f)) is met, this can be done with a for loop. Combined with the above recommendation, you might get something like this:
% Sort A in descending order instead, so first element in the array is now highest in value
[AVal, AIdx] = sort(A,1,'descend');
for f = 1:F
A_sorted_value = AVal(:,f); % Appears to be unused?
A_sorted_index = AIdx(:,f);
% Iterate through entire A_sorted_index array, eliminates rewriting of A_sorted_index
for i = 1:length(A_sorted_index)
v = A_sorted_index(i);
% Loop Exit Conditions
if (sum(OAK(1:t,f))<Yd(f)) || (sum(A_sorted_value(1:i)) <= 0)
break
end
% Your code below
%assign omega=1 if existing capacity is below desired production
if AK(v,f)+sum(OAK(1:t,f)) <= Yd(f)
omega(v,f)=1;
%update existing capacity
OAK(v,f)=omega(v,f).*AK(v,f);
%assign omega so that machine v can fill desired production
elseif AK(v,f)+sum(OAK(1:t,f))>Yd(f) && Yd(f)>sum(OAK(1:t,f))
omega(v,f) = (Yd(f)-sum(OAK(1:t,f)))/AK(v,f);
OAK(v,f) = omega(v,f).*AK(v,f);
end
end
end
Now, A_sorted_value appears to be unused. I can't assume that the value of elements in A are non-zero, positive numbers. So I had to include the check condition sum(A_sorted_value(1:i)) <= 0. But if you do have only non-zero, positive A elements, you can remove that snippet to speed the loop up even more.
I profiled this code against the old one, at T=200, F=100, ran in a loop 50 times.
new = 6.835 secs
old = 65.634 secs
Which is approximately a 90% time reduction.
A large amount of the time reduced came from removing these lines of codes, which were quite wasteful.
A_sorted_index(end)=[];
A_sorted_value(end)=[];
I'm sure there are still plenty more room for optimization, but this is good enough for me.
Do verify that the code I provided is correct.