时间戳

首页 > 解决方案 > 如何针对大量迭代优化这种随机熵计算？

arrays - 如何针对大量迭代优化这种随机熵计算？

问题描述

我的代码试图做什么：

我有一个包含分子排列的初始数组，每个分子都在一个单元格中，并且仅限于在二维数组中移动（上下左右）（数组大小：200x200）。在每一步，我都会随机抽取一个分子，并将其移动到一个随机的相邻单元格中。

从某个数字开始，每隔几次迭代，我计算这个网格的熵。网格被切割成 25x25 的小方块。然后我使用香农熵来计算系统的熵。

客观的：

使用 i5-6500，不使用 GPU，在适当的时间内进行 1e8+ 次迭代。

我的代码：

function advance_multi_lattice(grid)   #Find the next state of the system

    rnd=rand(1:count(!iszero,grid))        #Random number to be used for a random molecule.
    slots=find(!iszero,grid)    #Cells containing molecules.
    chosen_slot=find(!iszero,grid)[rnd]    #Random cell. May contain multiple molecules.
    dim=size(grid)[1]    #Need this for rnd=3,4 later.  
    grid[chosen_slot]-=1    #Remove the molecules from the cell
    rnd_arr=[1,2,3,4]     #Array to random from.


    while true


        rnd=rand(rnd_arr)   #Random number to see which side should the molecules go.

        if rnd==1    #Right for example.
            try     #In case moving right is impossible, ie: moving right gets the molecule out. Remove 1 from rnd_arr and repeat.
                grid[chosen_slot+1]+=1
                break
            catch
                filter!(e->e!=1,rnd_arr)
            end
        elseif rnd==2
            try   #Same
                grid[chosen_slot-1]+=1
                break
            catch
                filter!(e->e!=2,rnd_arr)
            end
        #Repeat for the other numbers : 3 and 4...
    return Grid
end

function S(P)   #Entropy, if no molecules then return 0.
    s=[]
    for k in P
        if k==0
            push!(s,0)
        else
            push!(s,-k*log(k))
        end

    end
    return s
end

function find_molecules(grid) #How many molecules in the array
    s=0

    for slot in grid
        s+=slot
    end
    return s
end

function entropy_scale(grid,total_molecules)    #Calculate the entropy of the grid.
    P_array=Array{Float64}([])
    for i=1:8
        for j=1:8
            push!(P_array,find_molecules(grid[(i-1)*25+1:i*25,(j-1)*25+1:j*25]))
        end
    end

    P_array=P_array./total_molecules

    return sum(S(P_array))
end

function entropy_evolution(grid,n)    #The loop function. Changes the grid and returns the entropy as a function of steps.
    t_arr=Array{Int64}([])
    S_arr=Array{Float64}([])
    p=Progress(Int(n))    #Progress bar, using ProgressMeter.
    total_molecules=find_molecules(grid)
    for k=1:1e3
        grid=advance_multi_lattice(grid)

        next!(p)
    end

    for k=1e3+1:n
        grid=advance_multi_lattice(grid)

        if  k%500==0  #Only record entropy every 500 steps
            push!(S_arr,entropy_scale(grid,totel_molecules))
        end
        next!(p)
    end

    return S_arr,grid

end

我的代码的结果：

对于 1e5 次迭代，我得到 43 秒。这意味着如果我想要一个有趣的结果（ 1e9+ ），我需要很多时间，最多 1 小时以上。更改熵计算阈值几乎不会影响性能，除非它非常小。

标签： arraysrandomjulia

解决方案

我假设您在 Julia 1.0 下工作（对于 Julia 0.6，需要进行一些小改动 - 我在代码中注明了这一点）。

为了提高性能，您应该保留分子向量 - 而不是网格（您不需要它，因为您允许分子占据相同的位置）。

我们将分子的位置编码为元组(x,y)。现在你需要一个随机移动一个分子的函数。以下是您如何实现它（我硬编码了边界，但当然您可以将它们更改为参数）：

function move_molecule((x,y)) # in Julia 0.6 it should be move_molecule(t)
    # and here in Julia 0.6 you should add: x, y = t
    if x == 1
        if y == 1
            ((1,2), (2,1))[rand(1:2)]
        elseif y == 200
            ((1,199), (2,200))[rand(1:2)]
        else
            ((2,y), (1,y-1), (1, y+1))[rand(1:3)]
        end
    elseif x == 200
        if y == 1
            ((200,2), (199,1))[rand(1:2)]
        elseif y == 200
            ((200,199), (199,200))[rand(1:2)]
        else
            ((200,y), (200,y-1), (200, y+1))[rand(1:3)]
        end
    else
        if y == 1
            ((x,2), (x-1,1), (x+1, 1))[rand(1:3)]
        elseif y == 200
            ((x,199), (x-1,200), (x+1, 200))[rand(1:3)]
        else
            ((x+1,y), (x-1,y), (x, y+1), (x,y-1))[rand(1:4)]
        end
    end
end

现在，一个将在一个步骤中移动给定数量的随机分子的函数steps是：

function go_sim!(molecules, steps)
    for k in 1:steps
        i = rand(axes(molecules, 1)) # in Julia 0.6 it should be: i = rand(1:length(molecules))
        @inbounds molecules[i] = move_molecule(molecules[i])
        if k % 500 == 0
        # here do entropy calculation
        end
    end
end

您没有提供完全可重现的示例，所以我在此停止 - 但使用此数据结构重写用于熵计算的其余代码应该很容易（实际上它可能更简单）。这是一个基准（性能不依赖于网格的大小，也不依赖于分子的数量，这是使用网格的代码的一个重要优势）：

julia> molecules = [(rand(1:200), rand(1:200)) for i in 1:1000];

julia> @time go_sim!(molecules, 1e9)
 66.212943 seconds (22.64 k allocations: 1.191 MiB)

你会1e9在大约一分钟内获得步骤（没有熵计算）。

良好性能所需的关键要素是什么：

  • 不要使用 try-catch 块，因为它们非常慢；
  • 尽量避免分配内存（即创建可变对象）；我的代码基本上没有分配 - 特别是这就是我到处使用元组的原因（move_molecule为简单起见，您可以在函数中使用矩阵，但性能会差 2 倍左右）

希望这可以帮助。

推荐阅读