python - 将无法向量化的 `for`循环转换为稀疏矩阵

Question

有 2 个盒子和一个小间隙，每秒允许 1 个粒子从一个盒子进入另一个盒子。粒子是否从 A 到 B，或从 B 到 A 取决于 Pa/Ptot 的比率(Pa:A 框中的粒子数，Ptot:两个框中的粒子总数)。

为了加快速度，I need to get rid of the for loops ，但是我找不到一种方法来向量化它们或将它们转换为表示我的 for 循环的稀疏矩阵:

What about for loops you can't vectorize? The ones where the result at iteration n depends on what you calculated in iteration n-1, n-2, etc. You can define a sparse matrix that represents your for loop and then do a sparse matrix solve.

但我不知道如何从中定义稀疏矩阵。模拟归结为计算:

哪里

当我试图表达我的问题时，这给我带来了麻烦 here 。 (注意:括号内的内容为 bool 运算)

问题:

• 我可以矢量化 for 循环吗？
• 如果没有，如何定义稀疏矩阵？
• (奖励问题)为什么 Python 中的执行时间(0.027 秒)比 Octave(0.75 秒)快x27？

_{注意:我在Python和Octave中实现了模拟，很快就会在Matlab上进行模拟，因此标签是正确的。}

<小时/>

Octave 代码

1; % starting with `function` causes errors

function arr = Px_simulation (Pa_init, Ptot, t_arr)
  t_size = size(t_arr);
  arr = zeros(t_size);   % fixed size array is better than arr = [] 
  rand_arr = rand(t_size);  % create all rand values at once
  _Pa = Pa_init;
  for _j=t_arr()
    if (rand_arr(_j) * Ptot > _Pa)
      _Pa += 1;
    else
      _Pa -= 1;
    endif
    arr(_j) = _Pa;
  endfor
endfunction


t = 1:10^5;

for _i=1:3
  Ptot = 100*10^_i;
  tic()
  Pa_simulation = Px_simulation(Ptot, Ptot, t);
  toc()
  subplot(2,2,_i);
  plot(t, Pa_simulation, "-2;simulation;")
  title(strcat("{P}_{a0}=", num2str(Ptot), ',P=', num2str(Ptot)))
endfor

<小时/>

Python

import numpy
import matplotlib.pyplot as plt
import timeit
import cpuinfo

from random import random

print('\nCPU: {}'.format(cpuinfo.get_cpu_info()['brand']))


PARTICLES_COUNT_LST = [1000, 10000, 100000]
DURATION = 10**5

t_vals = numpy.linspace(0, DURATION, DURATION)


def simulation(na_initial, ntotal, tvals):
    shape = numpy.shape(tvals)
    arr = numpy.zeros(shape)
    na_current = na_initial

    for i in range(len(tvals)):
        if random() > (na_current/ntotal):
            na_current += 1
        else:
            na_current -= 1
        arr[i] = na_current
    return arr


plot_lst = []
for i in PARTICLES_COUNT_LST:
    start_t = timeit.default_timer()
    n_a_simulation = simulation(na_initial=i, ntotal=i, tvals=t_vals)
    execution_time = (timeit.default_timer() - start_t)
    print('Execution time: {:.6}'.format(execution_time))
    plot_lst.append(n_a_simulation)


for i in range(len(PARTICLES_COUNT_LST)):
    plt.subplot('22{}'.format(i))
    plt.plot(t_vals, plot_lst[i], 'r')

    plt.grid(linestyle='dotted')
    plt.xlabel("time [s]")
    plt.ylabel("Particles in box A")

plt.show()

Answer 1

IIUC 您可以在 Octave 和 Numpy 中使用 cumsum():

Octave :

>> p = rand(1, 5);
>> r = rand(1, 5);
>> p
p =

   0.43804   0.37906   0.18445   0.88555   0.58913

>> r
r =

   0.70735   0.41619   0.37457   0.72841   0.27605

>> cumsum (2*(p<(r+0.03)) - 1)
ans =

   1   2   3   2   1

>> (2*(p<(r+0.03)) - 1)
ans =

   1   1   1  -1  -1

另请注意，以下函数将返回值 ([-1, 1]):