Subscript assignment dimension mismatch when calling function

I created a function that stochastically simulates cellular birth and death populations. When I copy and paste the function into a workspace, it runs fine. When I try to call the function in a livescript, I get the error "Subscript assignment dimension mismatch: line 191 pop_new(n,:) = population" This error only appears when the function is called, not when it is run separately. I have also checked to make sure the two vectors are the same size. They are both 1x7 doubles.

Here is the livescript code

mp_v = [1E-9, 1E-8, 1E-7];
mc_v = [1E-7, 1E-6, 1E-5];
pop_initial_v = [1E11, 0, 0, 0, 0, 0;
               1E10, 0, 0, 0, 0, 0;
               1E9, 0, 0, 0, 0, 0;
               1E8, 0, 0, 0, 0, 0];
mpc7_v = [1E-9, 1E-8, 1E-7];
f = zeros(108,10);
n = 1;
for i=1:length(mp_v)
  for j = 1:length(mc_v)
      for k = 1:length(pop_initial_v(1,:))
          for m = 1:length(mpc7_v)
              mp = mp_v(i);
              mc = mc_v(j);
              pop_initial = pop_initial_v(k,:);
              mpc7 = mpc7_v(m);
              [ end_pop ] = t790m(mp, mc, mpc7, pop_initial);
              f(n,:) = [n, mp, mc, pop_initial(1), mpc7, population(1), population(2), population(3), population(4), population(5)];
              n = n+1;
          end
      end
  end
end

and here is the function code

function [end_pop] = t790m(mp, mc, mpc7, pop_initial)
alpha = [.21, .21, .21, .21, .21, .21, .21]; %birth rate, for now assumed to be the same for all cells
beta = [.19, .19, .19, .19, .19, .19, .19]; %death rate before treatment
beta2 = [.23, .19, .19, .19, .19, .19, .19]; %death rate during gefitinib treatment
beta3 = [.23, .23, .21, .20, .20, .19, .19]; %death rate during PF00299804 treatment
tau_initial = NaN;
epsilon = .01;      %error tolerance
treatment_end = 1000;
initial_tauleap = true;
initial_hybrid = true;
initial_thresh = 1E5;
treatment1_start = 20;
treatment2_start = 300;
treatment2_end = 600;
%tau leaping method
t = 0;      %time vector
cells = pop_initial;
n = 1;      %counter
tau = tau_initial; 
singlessa = NaN;
while (sum(t)<=treatment2_end && sum(cells(end,:)>0))
  currtime = sum(t);
  population = cells(end,:);

    if currtime>=treatment1_start && currtime<treatment2_start
        gamma = beta2;  %death rate during treatment 1
    elseif currtime>=treatment2_start && currtime<=treatment_end
        gamma = beta3;   %death rate during treatment 2
    else
        gamma = beta;
    end

    events_s = [alpha(1)*(1-mp)*(1-mpc7)*population(1); %birth of s cells
                gamma(1)*population(1);         %death of s cell
                mp*alpha(1)*population(1);      %s cell mutates into r cell
                mpc7*alpha(1)*population(1)];   %s mutates into c797s
    state_change_s = [1, 0, 0, 0, 0, 0, 0;
                     -1, 0, 0, 0, 0, 0, 0;
                      0, 1, 0, 0, 0, 0, 0;  
                      0, 0, 0, 0, 0, 0, 1];%change in cellular populations due to each event
    events = zeros(25,1);
    events(1:4) = events_s;
    state_change = zeros(25,7);
    state_change(1:4,:) = state_change_s;
    if (initial_hybrid && n>2)
        if (population(1) < initial_thresh)
            rb = alpha;     rb(1) = 0;
            rd = gamma;     rd(1) = 0;
            events(1:2) = 0;
        end

        if (population(2) > 0)   %r cells have appeared
            events(5:8) = [alpha(2)*(1-mc)*(1-mpc7)*population(2);   %r replicates
                           gamma(2)*population(2);        %r dies
                           mc*alpha(2)*population(2);     %r mutates into rc1 cell
                           mpc7*alpha(2)*population(2)];                %r mutates into c797s
            state_change(5:8,:) = [0, 1, 0, 0, 0, 0, 0;
                                   0, -1, 0, 0, 0, 0, 0;
                                   0, 0, 1, 0, 0, 0, 0;
                                   0, 0, 0, 0, 0, 0, 1];
            if (population(2) > initial_thresh)
                rb = alpha;     rb(2) = 0;
                rd = gamma;     rd(2) = 0;
                events(5:6) = 0;
            end
        end

        if (population(3) > 0)   %rc1 cells have appeared
             events(9:12) = [alpha(3)*(1-mc)*(1-mpc7)*population(3);   %rc1 replicates
                             gamma(3)*population(3);      %rc1 dies
                             mc*alpha(3)*population(3);   %rc1 mutates into rc2 cell
                             mpc7*alpha(3)*population(3)];         %rc1 mutates into c797s
            state_change(9:12,:) = [0, 0, 1, 0, 0, 0, 0;
                                    0, 0, -1, 0, 0, 0, 0;
                                    0, 0, 0, 1, 0, 0, 0;
                                    0, 0, 0, 0, 0, 0, 1];
           if (population(3) > initial_thresh)
               rb = alpha;      rb(3) = 0;
               rd = gamma;      rd(3) = 0;
               events(9:10) = 0;
           end
        end

        if (population(4) > 0)   %rc2 cells have appeared
            events(13:16) = [alpha(4)*(1-mc)*(1-mpc7)*population(4);    %rc2 cell replicates
                             gamma(4)*population(4);        %rc2 cell dies
                             mc*alpha(4)*population(4);     %rc2 cell mutates into rc3 cell
                             mpc7*alpha(4)*population(4)];           %rc2 cell mutates into c797s
            state_change(13:16,:) = [0, 0, 0, 1, 0, 0, 0;
                                     0, 0, 0, -1, 0, 0, 0;
                                     0, 0, 0, 0, 1, 0, 0;
                                     0, 0, 0, 0, 0, 0, 1];
            if (population(4) > initial_thresh)
                rb = alpha;         rb(4) = 0;
                rd = gamma;         rd(4) = 0;
                events(13:14) = 0;
            end
        end

        if (population(5) > 0)  %rc3 cells have appeared
            events(17:20) = [alpha(5)*population(5);     %rc3 cells replicates
                             gamma(5)*population(5);     %rc3 cell dies
                             mc*alpha(5)*population(5);  %rc3 cell mutates into rc4 cell
                             mpc7*alpha(5)*population(5)];        %rc3 cell mutates into c797s
            state_change(17:20,:) = [0, 0, 0, 0, 1, 0, 0;
                                     0, 0, 0, 0, -1, 0, 0;
                                     0, 0, 0, 0, 0, 1, 0;
                                     0, 0, 0, 0, 0, 0, 1];
            if (population(5) > initial_thresh)
                rb = alpha;         rb(5) = 0;
                rd = gamma;         rd(5) = 0;
                events(17:18) = 0;
            end
        end

        if (population(6) > 0)   %rc4 cells have appeared
            events(21:23) = [alpha(6)*(1-mpc7)*population(6);   %rc4 cell replicates
                             gamma(6)*population(6);            %rc4 cell dies
                             mpc7*alpha(6)*population(6)];      %rc4 cell mutates into c797s cell
            state_change(21:23,:) = [0, 0, 0, 0, 0, 1, 0;
                                     0, 0, 0, 0, 0, -1, 0;
                                     0, 0, 0, 0, 0, 0, 1];
            if (population(6) > initial_thresh)
                rb = alpha;         rb(6) = 0;
                rd = gamma;         rd(6) = 0;
                events(21:22) = 0;
            end
        end

        if (population(7) > 0)  %c797s cells have appeared
            events(24:25) = [alpha(7)*population(7);    %c797s cell replicates
                             gamma(7)*population(7)];   %c797s cell divides
            state_change(24:25,:) = [0, 0, 0, 0, 0, 0, 1;
                                     0, 0, 0, 0, 0, 0, -1];
            if (population(7) > initial_thresh)
                rb = alpha;         rb(7) = 0;
                rd = gamma;         rd(7) = 0;
                events(24:25) = 0;
            end
        end
    end
    if (initial_tauleap && isnan(singlessa))
        if (isnan(tau))
    %partial derivatives from jacobian
            partial_derivative = [alpha(1)*(1-mp)*(1-mpc7), 0, 0, 0, 0, 0, 0; %S cells 
                                  gamma(1), 0, 0, 0, 0, 0, 0;      
                                  mp*alpha(1), 0, 0, 0, 0, 0, 0;
                                  mpc7*alpha(1), 0, 0, 0, 0, 0, 0;
                                  0, alpha(2)*(1-mc)*(1-mpc7), 0, 0, 0, 0, 0;  %r cells  
                                  0, gamma(2), 0, 0, 0, 0, 0;      
                                  0, mc*alpha(2), 0, 0, 0, 0, 0;
                                  0, mpc7*alpha(2), 0, 0, 0, 0, 0;
                                  0, 0, alpha(3)*(1-mc)*(1-mpc7), 0, 0, 0, 0;   %rc1 cells 
                                  0, 0, gamma(3), 0, 0, 0, 0;       
                                  0, 0, mc*alpha(3), 0, 0, 0, 0;
                                  0, 0, mpc7*alpha(3), 0, 0, 0, 0;
                                  0, 0, 0, alpha(4)*(1-mc)*(1-mpc7), 0, 0, 0;    %rc2 cells   
                                  0, 0, 0, gamma(4), 0, 0, 0;         
                                  0, 0, 0, alpha(4)*mc, 0, 0, 0;
                                  0, 0, 0, alpha(4)*mpc7, 0, 0, 0;
                                  0, 0, 0, 0, alpha(5)*(1-mc)*(1-mpc7), 0, 0;  %rc3 cells      
                                  0, 0, 0, 0, gamma(5), 0, 0; 
                                  0, 0, 0, 0, alpha(5)*mc, 0, 0;
                                  0, 0, 0, 0, alpha(5)*mpc7, 0, 0;
                                  0, 0, 0, 0, 0, alpha(6)*(1-mpc7), 0;    %rc4 cells      
                                  0, 0, 0, 0, 0, gamma(6), 0;
                                  0, 0, 0, 0, 0, alpha(6)*mpc7, 0;
                                  0, 0, 0, 0, 0, 0, alpha(7);       %c797s cells
                                  0, 0, 0, 0, 0, 0, gamma(7)];
            total_rate = sum(events);
            fjj = partial_derivative * state_change';
            mu = fjj*events;    %mean
            std = (fjj.^2)*events;     %standard deviation
            tau = min(min(epsilon*total_rate./abs(mu), (epsilon^2)*total_rate^2./std));
        end
        if (tau < 1/total_rate)
            singlessa = 10;
            tau = tau_initial;
            continue
        else
            pop_new = zeros(800,7);
            pop_new(:,:) = 1E14;
            pop_new(n,:) = population;
            pop_new(n,:) = pop_new(n,:) + (state_change'*poissrnd(events*tau))'; %state_change'*poissrnd(events*tau)new population sizes 
            pop_new(n,:) = round(pop_new(n,:));
            if sum(pop_new(n)<0)>0
                tau = tau/2;
                continue
            else
                time_elapsed = tau;
                t(n) = time_elapsed;
                cells(n,:) = pop_new(n,:);
                tau = tau_initial;
                n = n+1;
            end
        end
    else
        r1 = rand;
        time_elapsed = -log(r1)./total_rate_s;
        if (total_rate > 0)
            t(n) = time_elapsed;
            idx = randsample(length(events_s),1,true,events_s);
            cells(n,:) = round(population + state_change_s(idx,:));
            n = n+1;
        end

        if (initial_tauleap)
            singlessa = singlessa - 1;
            if (singlessa < 1)
                singlessa = NaN;
            end
        end
    end
    if (currtime + time_elapsed > treatment2_end && time_elapsed > .1*treatment2_end)
        tspancalc = treatment2_end - sum(t(1:end-1));
        net = alpha(1) - gamma(1);
        cells(end,1) = cells(end-1,1)*exp(net*tspancalc);
        net = alpha(2) - gamma(2);
        cells(end,2) = cells(end-1,2)*exp(net*tspancalc);
        net = alpha(3)-gamma(3);
        cells(end,3) = cells(end-1,3)*exp(net*tspancalc);
        net = alpha(4)-gamma(4);
        cells(end,4) = cells(end-1,4)*exp(net*tspancalc);
        net = alpha(5)-gamma(5);
        cells(end,5) = cells(end-1,5)*exp(net*tspancalc);
        net = alpha(6)-gamma(6);
        cells(end,6) = cells(end-1,6)*exp(net*tspancalc);
        net = alpha(7) - gamma(7);
        cells(end,7) = cells(end-1,7)*exp(net*tspancalc);
        t(end) = tspancalc;
    end
    if (initial_hybrid && n>2)
        time_elapsed = t(end);
        if (population(1) > initial_thresh)
            net = alpha(1)-gamma(1);
            cells(end,1) = cells(end,1)+cells(end-1,1)*(exp(net*time_elapsed)-1);
        end

        if (population(2) > initial_thresh)
            net = alpha(2)-gamma(2);
            cells(end,2) = cells(end,2)+cells(end-1,2)*(exp(net*time_elapsed)-1);
        end

        if (population(3) > initial_thresh)
            net = alpha(3)-gamma(3);
            cells(end,3) = cells(end,3)+cells(end-1,3)*(exp(net*time_elapsed)-1);
        end

        if (population(4) > initial_thresh)
            net = alpha(4)-gamma(4);
            cells(end,4) = cells(end,4)+cells(end-1,4)*(exp(net*time_elapsed)-1);
        end

        if (population(5) > initial_thresh)
            net = alpha(5)-gamma(5);
            cells(end,5) = cells(end,5)+cells(end-1,5)*(exp(net*time_elapsed)-1);
        end

        if (population(6) > initial_thresh)
            net = alpha(6)-gamma(6);
            cells(end,6) = cells(end,6)+cells(end-1,6)*(exp(net*time_elapsed)-1);
        end

        if (population(6) > initial_thresh)
            net = alpha(7)-gamma(7);
            cells(end,7) = cells(end,7)+cells(end-1,7)*(exp(net*time_elapsed)-1);
        end
    end
end

t = cumsum(t'); t = [0; t]; cells = [pop_initial; cells]; end_pop = population;

%s_cells = population(1); %r_cells = population(2); %rc1_cells = population(3); %rc2_cells = population(4); %c797s = population(7) end

Answer 1

[James Tursa]

Type the following at the MATLAB prompt

dbstop if error

Then run your code. When the error occurs, the code will pause at the offending line with all current variables still intact. Examine them to see what the dimensions really are, and then you can backtrack to figure out why the dimensions are not what you expected.

[Kate Yorko]

That worked! thank you