help with parallel computing application

[dukandender]

hello, i am performing research in which i need to scale up my slow function evaluation massively. i am performing optimization with BlackBoxOptim.jl and would like to parallelize this in the following way with Dagger and DTables:

1. All workers learn all the functions required to eventually compute loss(), which is the function to be optimized. I plan to do this with @Everywhere.
2. Asynchronously have workers individually search the solution space with probablistic descent for 20 seconds with an initial guess x0. The two outputs (namely the optimizable parameters and the loss/fitness) will be saved for comparison.
3. As each worker completes its 20 second run, the best loss/fitness value will be compared to the absolute best loss/fitness value amongst all workers. I want to use a shared storage (potentially DTables?) so that each worker can access the current best loss/fitness value for comparison, and update it if necessary. I need this to be locked/mutexed so that comparison happens one-by-one, and the absolute lowest parameter-loss pair is saved properly.
4. As the shared table is altered by the workers optimizing their functions, each new 20 second run starts by updating x0 to be the current best parameter set (defined by the best/lowest loss/fitness value).
5. In case the program crashes and I'd like a backup, I would like to save the current best parameter-loss pair to a txt file with writedlm() every four cycles. This would only happen on one thread to prevent IO issues.
6. I want to scale this up to many, many threads, as the function evaluation is quite slow and I would like the optimization to benefit from massively searching the solution space. Ideally, I would want to be able to parallelize this on a GPU, as I have access to a HPC with A100s.

For my current setup, the probablistic descent works very well and seemingly always decreases on a 20second cycle. However, this drop is extremely minimal, so by paralellizing, hopefully this optimization will happen much, much quicker. To restate my aim, I would like to parallelize my optimizer by having workers search the space concurrently and properly handle saving the best pairs, when they come up.

I am new to Dagger and DTables and Distributed, so any help would be extremely helpful. If anyone is able to help me conceptualize this in a MWE, I will be able to probably adapt this for my application. For now, here is how I want to call the optimizations, but if a function is not the best way to allow for parallel computing, please let me know. Thank you for reading this and let me know if anything needs to be clarified.
julia
function optimize_worker(initial_params, duration)
res = bboptimize(
loss, initial_params;
NumDimensions=length(initial_params),
MaxTime=duration,
SearchRange=(-17, 17),
TraceMode=:silent,
PopulationSize=5000,
Method=:probabilistic_descent,
lambda=100,
)
best_params = best_candidate(res)
best_fitness = best_fitness(res)
return (best_params=best_params, best_fitness=best_fitness)
end

[jpsamaroo]

Hi @dukandender ! In order to help you determine how best to use Dagger to parallelize your evaluations, it would be helpful to know what kind of function you're evaluating, especially if you want to run it on GPUs.

I would personally recommend a simpler approach than using DTables, which would be just using Dagger.@spawn myfunc(...) together with using a ReentrantLock-protected struct or DataFrame. This will be simple to write and maintain, and since your function takes a long time to execute, you shouldn't spend any significant amount of time with threads contending on the lock (just keep the logic within the locked region as simple and fast as possible).

As an example of an inner loop:

N_evals = 100

lck = Threads.ReentrantLock()
current_best_params = ...
current_best_fitness = ...
function update_and_get_best!(result)
  lock(lck) do
    if best_fitness(result) < current_best_fitness
      current_best_fitness = best_fitness(result)
      current_best_params = best_candidate(result)
    end
    # Now could be a good time to checkpoint your current best result!
    # Consider incrementing a `Ref{Int}` counter here to determine when to checkpoint
    return (current_best_fitness, copy(current_best_params))
  end
end

function run_sim!()
  local_best_fitness = ...
  local_best_candidate = ...
  for eval_idx in 1:N_evals
    result = Dagger.@spawn myfunc(local_best_candidate)
    wait(Dagger.@spawn scope=Dagger.scope(worker=myid()) update_and_get_best!(result))
  end
end

@sync for _ in 1:Threads.nthreads()
  Threads.@spawn run_sim!()
end
# Final result is in `current_best_*`

[dukandender]

Hey Julian! Thanks for the reply. I've been trying lots of things since I made the post initially. Before I implement your approach, I would like your advice on a few things, as I think my questions have shifted slightly. Here is the way I’ve been doing things currently:

while true
    local pd = vec(readdlm("out.txt")) # pd are my params
    local curloss = consolidatedLoss(pd) # this is my current loss
    t1 = time()
    tasks = [Distributed.@spawn opt(consolidatedLoss, pd) for _ in 1:(nproc()/2)]
    #note in the line above, opt() is a function I've defined somewhere else that runs Optim's PSO with my set criteria (e.g., iterations, bounds, etc.)
    #also note that i have to do nproc()/2, not nproc(), to get accurate finishing times … do you know why that is? (e.g., if i set time limit to 7sec on my Optim.jl PSO optimizer, it will come out as 14sec if I do nproc() and 7sec if I do nproc()/2).
    results = map(Distributed.fetch, tasks)
    indecks = argmin([Optim.minimum(results[i]) for i in 1:length(results)])
    
    local pd = Optim.minimizer(results[indecks])
    best = Optim.minimum(results[indecks])
    curloss > best ? writedlm("out.txt", pd) : nothing
end

However, there are a few issues with this. First of all, loading the CPUs takes a while. For 400 procs, it took about 30 minutes. As my compute cluster uses LSF/bsub, I do:
JULIA_DEPOT_PATH="/home/research/me/.julia" JULIA_WORKER_TIMEOUT='1000' julia --project="/storage1/me/Projects/myproj" --procs 400 src/main.jl
Note that I need to set the timeout to 1000, or else it breaks. Here are my questions:

1. Comparing my approach and your approach, is it worth it to rework my code at this stage? What are the benefits of using Dagger/your approach for my application, outside of maintaining the code?
2. It’s seemed that I need to do nprocs()/2 (see commented note in the code). Why do I need to do this to parallelize as expected? Am I making a mistake somewhere and I don't need to do this?
3. In your approach, I don't have the conceptual understanding of why it makes sense to do Dagger.@Spawn within the loop for each thread. I may not fully understand the difference between a thread and a proc. What specifically does Dagger.@Spawn allow for, if it is being run on one thread? I guess I've always thought about it like multiple "computers" performing the same task at the same time, and then congregating after the task is over. When you iterate over many threads, I assume you're starting the task on each of those computers, so what exactly is Dagger.@Spawn doing, at a high level?

And probably the most important question I have:
5. Eventually, I would like to use PSOGPU.jl to optimize my function, as there are 168 parameters (soon, there will be a LOT more, like around 5k, after I show that it optimizes well for my current case) and the loss function takes around 20ms on benchmark. As it stands, waiting 30 min to start up my parallelizing, cutting my usable procs in half, and the CPU usage limits of my cluster all have dissuaded me from optimizing with CPU parallelization. If I were to do this, I would need to make my loss function GPU callable, right? How would I go about checking this? I've already started by implementing StaticArrays, but I still get a few errors when running PSOGPU ("InvalidIRError" and "invalid LLVM IR
Reason: unsupported call to an unknown function (call to julia.new_gc_frame)
Reason: unsupported call to an unknown function (call to julia.push_gc_frame)
Reason: unsupported call to an unknown function (call to julia.get_gc_frame_slot)", which I am assuming are due to other functions that I am calling within my master loss function that have inputs not prepared for GPU computation). I can make a new discussion thread on the PSOGPU Github. Or, if you think it would be more efficient to optimize with the CPU, then I can work on making my current process more efficient. Thank you!

[pszufe]

@dukandender I do optimization and got thinking about your problem statement. Typically in optimization you have a global state and need to iterate over many small jobs that rely on that global state, Than the state needs to be updated and used in subsequent runs of those jobs.
Additionally computations performed on the global state are very simple (contrary to the long running jobs).

In such optimization setting another approach would be to use green threads to dispatch remote jobs, but manipulate the global state in a single thread. The code could look more or less like this:

using Distributed, DataFrames,Dates
addprocs(4)
@everywhere using Random, Distributions
@everywhere Random.seed!(myid())

# of course you add to the global state much more info 
@kwdef mutable struct OptState
    log::DataFrame=DataFrame(step=Int[], wid=Int[], starttime=DateTime[], finish=DateTime[], x=Float64[], val=Float64[])
    bestix::Int=0
    bestixwid::Int=0
    bestixstep::Int=0
    freeWorkers::Set{Int}=Set(workers())
end

@everywhere f(x) = exp(sin(x))*2x^3+exp(cos(x))*2x^2-12x+4+mean(rand(Normal(0,1), rand(1:1000))) # some random function

@everywhere function optimize(f)
    bst = argmax(f, rand(100_000).*100 .- 50)
    bst, f(bst)
end

os = OptState()
asyncmap(1:100;ntasks=nworkers()) do step
    # this gets executed always in the main process
    # and single threaded so we can use the global state
    starttime = now()
    wid = pop!(os.freeWorkers)
    # you can desing optimize in a way to pass any parameters precomputed from the global state `os`
    x, val = fetch(@spawnat wid optimize(f)) # this is blocking io so this will yield
    # the lines below are no-blocking so they will not yield - will be executed in a single continues block so no race issues
    push!(os.freeWorkers, wid)
    push!(os.log, (;step, wid, starttime, finish=now(), x, val))
    if os.bestix == 0 || val > os.log.val[os.bestix]
        os.bestix = nrow(os.log)
        os.bestixwid = wid
        os.bestixstep = step
    end
end

Does this template look like useful for your purpose?