Fit quality parameter in a double loop

I have a physical process (amplification in gain medium) described by the equation for the parameter "amp_I0" in the code below. The free parameters in the process are the variables g0 and alpha, and the goal is to calculate the residual sum of squares (fitquality parameter in the code) between the reference and amplified profiles. Everything is working fine when I manually change g0 and alpha (the commented section in the code), but when I want to change the free parameters in loops and get the fitquality parameter for each g0 and alpha, I get an array mostly populated with zeros (which does not make sense in terms of the physcis). Can someone please fix my mistake in the double for loop? The strange thing is the plot looks correct.

clear
close all

t = [-100:0.1:100];
I = exp(-(t/10).^2) + 0.5*exp(-((t-20)/10).^2);
I = I/max(I);
I = I';

ref_I = exp(-(t/4.5).^2) + 0.06*exp(-((t-20)/6).^2);
ref_I = ref_I/max(ref_I);
ref_I = ref_I';

lambda = [190:0.01:197];
S = exp(-((lambda-193.5)/0.25).^2);
S = S/max(S);
S = S';

L = 54;
g0 = 1*0.16;
alpha = 10*0.008;

prod = S*I';

% amp_I0 = g0*prod.*(exp((g0*prod-alpha)*L)-1)./(g0*prod-alpha);
% amp_I = sum(amp_I0,1)/max(sum(amp_I0,1));
% amp_I = amp_I';
% amp_I = amp_I/max(amp_I);
% fitquality = trapz((amp_I - ref_I).^2)

for k = 1:11
    alpha(k) = 0.008 + 0.008*(k-1);
    for j = 1:21
        g0(j) = 0.15 + 0.001*(j-1);
        amp_I0 = g0(j)*prod.*(exp((g0(j)*prod-alpha(k))*L)-1)./(g0(j)*prod-alpha(k));
        amp_I(:,j*k) = sum(amp_I0,1)/max(sum(amp_I0,1));
        amp_I(:,j*k) = amp_I(:,j*k)';
        amp_I(:,j*k) = amp_I(:,j*k)/max(amp_I(:,j*k));
        fitquality(k,j*k) = trapz((amp_I(:,j*k) - ref_I).^2);
        % plot(t,amp_I,'r');
        % drawnow()
    end
end

figure(1)
plot(t,I,'b')
hold on
plot(t,ref_I,'r')
plot(t,amp_I,'g')

figure(2)
plot(lambda,S,'b')

Answer 1

[Voss]

The problem is using j*k as the index where you store the results. Consider when k=1, j*k is 1,2,3,4,5,6,...,21. And when k=2, j*k is 2,4,6,8,10,12,...,42. Notice the list for k=2 includes some values of j*k that were already used when k was 1. Thus, some results from k=1 are overwritten when k=2 (and some of those will be overwritten again when k is 3, and so on).

[You already know that 11*21 is the number of (j,k) pairs you have, so that's how many columns amp_I should have, and that's how many iterations you have in the code, which is good. But since some columns are being written to more than once as shown in the previous paragraph, then that means some are not being written to at all, and in those columns MATLAB uses zeros as placeholders, which is why you see lots of columns of all zeros.]

All that is to say that j*k is not the correct expression for the column of amp_I where the results for (j,k) should be stored. Proper expressions include (k-1)*21+j and (j-1)*11+k, depending on the order you want.

[Voss]

And of course it's better not to hard-code numbers like 11 and 21 into your expressions, but instead store them as variables, say nk=11; nj=21; and use those variables throughout your code, e.g., (k-1)*nj+j or (j-1)*nk+k.

Also, Pre-allocating amp_I before your loops is best practice. amp_I=zeros(numel(prod),nk*nj); would do it.

But you shouldn't name a variable "prod" because that name is already in use by a pretty important builtin function.

[Aram Zeytunyan]

This was exactly what I needed, thank you! Pre-allocating all the arrays inside the loop also saved me 3 seconds per each run.

[Voss]

You're welcome!