Rewrite Matlab ODE into Python - IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

How can I rewrite the following ODE calculation code into Python code?

countries = 150;
verbose = false;
export_raw_data = true;
if verbose
    countries=1;
end
if export_raw_data
    raw_features = zeros(4*countries, 501);
    raw_features(1,:) = 1:501;
end
features = zeros(countries,8);
rng('default')  % For reproducibility
betas = normrnd(0.0005,0.0001,[1,countries]);
for i=1:countries
    res = simVirusSpreading(betas(i), false);
    times = res(:,1);
    % Add randomness to the observations
    rand_res = zeros(size(res,1),5);
    rand_res(:,1) = times;
    for j=2:5
        rand_res(:,j) = arrayfun(@(x) normrnd(0,0.025*x)+x, res(:,j));
    end
    if verbose
        fprintf("Beta = %.3f%%\n",betas(1)*100);
        plot(times,[res(:,2) rand_res(:,2), res(:,3) rand_res(:,3),res(:,4) rand_res(:,4),res(:,5) rand_res(:,5)]);        
    end
    if export_raw_data
        days = [0:500];
        for var = 2:5
            F = griddedInterpolant(rand_res(:,1),rand_res(:,var));
            raw_features((i-1)*4+var,:)=F(days);
        end
    end
    ti = times(times<=50);
    n_50 = rand_res(length(ti),3);
    d_50 = rand_res(length(ti),3) - rand_res(length(ti)-1,3);
    ti = times(times<=150);
    n_150 = rand_res(length(ti),3);
    d_150 = rand_res(length(ti),3) - rand_res(length(ti)-1,3);
    ti = times(times<=300);
    n_300 = rand_res(length(ti),3);
    d_300 = rand_res(length(ti),3) - rand_res(length(ti)-1,3);
    [M,ii] = max(rand_res(:,3));
    t_peak = times(ii);
    if verbose
        fprintf("n_50=%.4f, n_150=%.4f, n_300=%.4f\nd_50=%.4f, d_150=%.4f, d_300=%.4f\nt_peak=%.4f\nx_500=%.4f\n",n_50, n_150, n_300, d_50, d_150, d_300, t_peak, rand_res(end,5) );
    end
    features(i,:)=[n_50, n_150, n_300, d_50, d_150, d_300, t_peak, rand_res(end,5)];
end
writematrix(features,'../2020 DL data/epidemic_process.csv');
if export_raw_data
    writematrix(raw_features,'../2020 DL data/epidemic_process_raw_data.csv');
end

My approach to write the Matlab code above into Python

from scipy.interpolate import RegularGridInterpolator

countries = 150
verbose = False
export_raw_data = True

if verbose:
    countries = 1
    
if export_raw_data:
    raw_features = np.zeros([4 * countries, 501])
    raw_features[0,:] = np.arange(1, 502)

features = np.zeros([countries, 8])
betas = np.random.default_rng().normal(0.0005, 0.0001, [1, countries])

for i in range(0, countries):
    res = simVirusSpreading(betas[0][i], False)
    times = res[:,0]
    rand_res = np.zeros([res.shape[0], 5])
    print(rand_res.shape)
    rand_res[0] = times[:,0]
   
    for j in range(1, 5):
        rand_res[:,j] = (lambda x: normrnd(0,0.025*x)+x, res[:,j]);
        #[tSol, ySol] = ode45(@(t,y) virusSpreading(t, y, beta, mu, nu), t_span, y0);
        #sol = solve_ivp(lambda t, y: virusSpreading(t, y, beta, mu, nu), t_span, y0);
        ##rand_res(:,j) = arrayfun(@(x) normrnd(0,0.025*x)+x, res(:,j));
        
    
    if verbose:
        ##fprintf("Beta = %.3f%%\n",betas(1)*100);
        plt.plot(times,[res[:,1], rand_res[:,1], res[:,2], rand_res[:,2],res[:,3], rand_res[:,3],res[:,4], rand_res[:,4]]); 

    if export_raw_data:
        days = np.arange(0, 501)
        
        for var in range(1,5):
            F = RegularGridInterpolator(rand_res[:,0], rand_res[:,var])
            raw_features[(i-1)*4+var,:] = F[days]
            
    ti = times[times<=50]
    n_50 = rand_res[(ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1]), 2]
#   d_50 = rand_res(length(ti),3) - rand_res(length(ti)-1,3);

    d_50 = rand_res[(ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1]), 2] - rand_res[ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1] -1, 3]
    ti = times[times<=150]
    n_150 = rand_res[(ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1]), 2]
    d_150 = rand_res[(ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1]), 2] - rand_res[(ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1])-1, 3]
    ti = times[times<=300]
    n_300 = rand_res[(ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1]), 2]
    d_300 = rand_res[ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1], 2] - rand_res[(ti.shape[0] if ti.shape[0] >= ti.shape[1] else ti.shape[1]) -1, 3]

    M, ii = rand_res[:, 2].max(axis=0)
    t_peak = times[ii]
    
    if verbose:                                                                                                 
        ##fprintf("n_50=%.4f, n_150=%.4f, n_300=%.4f\nd_50=%.4f, d_150=%.4f, d_300=%.4f\nt_peak=%.4f\nx_500=%.4f\n",n_50, n_150, n_300, d_50, d_150, d_300, t_peak, rand_res(end,5) );    
        continue
    features[i,:] = [n_50, n_150, n_300, d_50, d_150, d_300, t_peak, rand_res[end,4]]

np.savetxt("../epidemic_process.csv", features, delimiter=",")

if export_raw_data:
    np.savetxt("../epidemic_process.csv", raw_features, delimiter=",")
    
    

I always receive in Python the following code line:

IndexError                                Traceback (most recent call last)
<ipython-input-145-fced60228e63> in <module>
     73 for i in range(0, countries):
     74     res = simVirusSpreading(betas[0][i], False)
---> 75     times = res[:,0]
     76     rand_res = np.zeros([res.shape[0], 5])
     77     print(rand_res.shape)

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

If you need further information about the "simVirusSpreading" function, please let me know and I will provide the code for it as well.

Answer 1

[Jan]

The error message tells, that res replied by simVirusSpreading() is a vector, but you try to access it as a matrix with 2 indices. I have too few experiences in Python and this is a Matlab forum, but with pure guessing I'd try this:

times = res[:]