Parallel the particle swarm optimization in Julia

Question

I wrote the particle swarm optimization algorithm in Julia, and I tried to use Threads to accelerate the computation for each population. Here is the code

function uniform(dim::Int, lb::Array{Float64, 1}, ub::Array{Float64, 1})
    arr = rand(Float64, dim)
    @inbounds for i in 1:dim; arr[i] = arr[i] * (ub[i] - lb[i]) + lb[i] end
    return arr
end

function initialize_particles(problem, population)
    dim = problem.dim
    lb = problem.lb
    ub = problem.ub
    cost_func = problem.cost_func

    gbest_position = uniform(dim, lb, ub)
    gbest = Gbest(gbest_position, cost_func(gbest_position))

    particles = []
    for i in 1:population
        position = uniform(dim, lb, ub)
        velocity = zeros(dim)
        cost = cost_func(position)
        best_position = copy(position)
        best_cost = copy(cost)
        push!(particles, Particle(position, velocity, cost, best_position, best_cost))

        if best_cost < gbest.cost
            gbest.position = copy(best_position)
            gbest.cost = copy(best_cost)
        end
    end
    return gbest, particles
end

function PSO(problem; max_iter=100,population=100,c1=1.4962,c2=1.4962,w=0.7298,wdamp=1.0)
    dim = problem.dim
    lb = problem.lb
    ub = problem.ub
    cost_func = problem.cost_func

    gbest, particles = initialize_particles(problem, population)

    # main loop
    for iter in 1:max_iter
        @threads for i in 1:population
            particles[i].velocity .= w .* particles[i].velocity .+
                c1 .* rand(dim) .* (particles[i].best_position .- particles[i].position) .+
                c2 .* rand(dim) .* (gbest.position .- particles[i].position)

            particles[i].position .= particles[i].position .+ particles[i].velocity
            particles[i].position .= max.(particles[i].position, lb)
            particles[i].position .= min.(particles[i].position, ub)

            particles[i].cost = cost_func(particles[i].position)

            if particles[i].cost < particles[i].best_cost
                particles[i].best_position = copy(particles[i].position)
                particles[i].best_cost = copy(particles[i].cost)
        
                if particles[i].best_cost < gbest.cost
                    gbest.position = copy(particles[i].best_position)
                    gbest.cost = copy(particles[i].best_cost)
                end
            end
        end
        w = w * wdamp
        if iter % 50 == 1
            println("Iteration " * string(iter) * ": Best Cost = " * string(gbest.cost))
            println("Best Position = " * string(gbest.position))
            println()
        end
    end
    gbest, particles
end

Obviously, it is not thread-safe in main loop of PSO. This is because the position of gbest needs to be modified when computing. I tried to use Atomic, but the position of each particle is an array. So this method is not suitable. The Lock needs initialization, so it also won't work. I also tried to use Floop, but it went wrong.

So is there a way to make sure the data is locked when doing parallel computation?

Thank you for your time!

Wang

Andrej Oskin Andrej Oskin · Accepted Answer · 2020-11-10T12:58:15

There is no need to lock variables. When you are doing parallel calculations you can always think in the terms of Split-Combine strategy (sorry just invent this term in analogy with DataFrames calculation, you may also call it Map-Reduce approach). The idea is that you split calculations around different threads and do them independently and on the final stage you combine results together on a single thread. So your code can look like this (there can be syntax errors, since I can't run the code without Particle and other definitions, but I hope it is enough to give the idea)

function minby(itr; by=identity, init=nothing)
     init = isnothing(init) ? pop!(itr) : init
     mapreduce(x->(by(x)=>x), (x,y)->(first(x)>first(y)) ? y : x, itr; init=by(init)=>init).second
end

function PSO(problem; max_iter=100,population=100,c1=1.4962,c2=1.4962,w=0.7298,wdamp=1.0)
    dim = problem.dim
    lb = problem.lb
    ub = problem.ub
    cost_func = problem.cost_func

    gbest, particles = initialize_particles(problem, population)

    # main loop
    for iter in 1:max_iter
        gbests = fill((gbest.cost, 0), Threads.nthreads())
        @threads for i in 1:population
            particles[i].velocity .= w .* particles[i].velocity .+
                c1 .* rand(dim) .* (particles[i].best_position .- particles[i].position) .+
                c2 .* rand(dim) .* (gbest.position .- particles[i].position)

            particles[i].position .= particles[i].position .+ particles[i].velocity
            particles[i].position .= max.(particles[i].position, lb)
            particles[i].position .= min.(particles[i].position, ub)

            particles[i].cost = cost_func(particles[i].position)

            if particles[i].cost < particles[i].best_cost
                particles[i].best_position = copy(particles[i].position)
                particles[i].best_cost = copy(particles[i].cost)
        
                if particles[i].best_cost < gbests[Threads.threadid()][1]
                    gbests[Threads.threadid()] = (particles[i].best_cost, i)
                end
            end
        end
        gbest1 = minby(gbests, by = x -> x[1])
        if gbest1[2] != 0
            idx = gbest1[2]
            gbest.position = copy(particles[idx].best_position)
            gbest.cost = copy(particles[idx].best_cost)
        end
        w = w * wdamp
        if iter % 50 == 1
            println("Iteration " * string(iter) * ": Best Cost = " * string(gbest.cost))
            println("Best Position = " * string(gbest.position))
            println()
        end
    end
    gbest, particles
end

By the way, you may find that package UnPack.jl is rather convinient. Instead of manual assignments, you can do just

using UnPack
function PSO(problem; max_iter=100,population=100,c1=1.4962,c2=1.4962,w=0.7298,wdamp=1.0)
    @unpack dim, lb, ub, cost_func = problem
    ....

Parallel the particle swarm optimization in Julia

1 Answers