Parallell map for ODE problems

[halvorak]

I want to solve an initial value problem (IVP) many times, for different parameters (x0 is constant). Now, I do this by the map function, see the example code below (CasADi version 3.5.5). I've two questions

1. How do I specify that e.g. openmp should use more threads? I can't find any documentation about this. I can specify e.g. F2.map(N, "thread", 4) and then use 4 threads, but if I use F2.map(N, "openmp", "thread", 4) does not work. If I specify F2.map(N, "openmp") it does not seem that any parallell computing is done. Basically - how do I pass extra info to the map construct?
2. Is it possible to solve IVP in parallell directly by using CVODES? In Casadis documentation for cvodes it mentions that "CVODES can be used both on serial and parallel computers. The main difference is in the NVECTOR module of vector kernels. The desired version is obtained when compiling the example files by linking the appropriate library of NVECTOR kernels. In the parallel versions,communication between processors is done with MPI, with OpenMP, or with Pthreads." I cannot find out how to evaluate IVPs in parallel without making a map construct though.

`if name == "main":
import numpy as np
import casadi as ca

# Declare variables
x  = ca.MX.sym('x', 3)  # state
p  = ca.MX.sym('p', 1)  # parameters

# ODE right hand side function and integrator
xdot = ca.vertcat((1 - x[1]*x[1])*x[0] - x[1] + p, \
              x[0], \
              x[0]*x[0] + x[1]*x[1] + p*p)
ode = {'x':x, 'p':p, 'ode':xdot}
F = ca.integrator("F", "cvodes", ode)

#initial conditions
x0 = np.array([0., 1., 0])

#functions which takes p as input. To be used in map
F_out = F(x0 = x0, p = p)["xf"]
F2 = ca.Function("F2", [p], [F_out])

#define parameter samples
N = int(1e4)
par_samples = np.random.normal(.2, .2*.01, size = N).reshape(1,-1)

#define map functions
F_map = F2.map(N)
F_map_t = F2.map(N, "thread", 4)
F_map_openmp = F2.map(N, "openmp", 4)
F_map_serial = F2.map(N, "serial")
F_map_unroll = F2.map(N, "unroll")

#run time
print(f"F_map, N = {N}")
%timeit sol = F_map(par_samples)
print(f"F_map_t")
%timeit sol = F_map_t(par_samples)
print(f"F_map_openmp")
%timeit sol = F_map_openmp(par_samples)
print(f"F_map_serial")
%timeit sol = F_map_serial(par_samples)
print(f"F_map_unroll")
%timeit sol = F_map_unroll(par_samples)

`

[brunoml5]

Hi! I remember some discussions on the old forum about openmp. I don't remember very well the details, but it appers that openmp is not working properly for a while. Only thread works well. I'm not a developer, so I don't know the inner workings of this feature, and the documentation is not very clear. Hopefully they answer so that we can learn.

[jaeandersson]

1. To use the OpenMP support, you CasADi to have been compiled with OpenMP support. The number of threads it will use will depend on how many threads are available in the system. If openmp isn't working properly, as @brunoml5 suggests, maybe try to use thread instead.
2. The parallelization support in Sundials/CVODES you mention is not relevant in this context. We only support serial NVECTOR modules in our Sundials interface. CasADi was developed with focus on optimization and then it's usually a better idea to parallelize at the optimizer level and not at the integrator level. Typically, you would do this by implement an algorithm like collocation or multiple shooting, which can easily be parallelized at the optimizer level and also have other benefits in terms of faster convergence.

[halvorak]

Thanks for the reply, and I fully agree on 2). I was just curious since NVECTOR was mentioned in the documentation.

I used CasADi version 3.5.5 in my original question. Now, I switched to version 3.6.2 (Python 3.10.11 on Windows, packaged by Anaconda). When I updated, thread does not work any longer (so now, all the map approaches in my original question evaluates the ODE in serial). Is there a guide on how to compile CasADi with openmp somewhere?

PS: on CasADi 3.6.2, I get the following warning about compilation of openmp. I did not get this warning on version 3.5.5 (so it is good that this warning is here in 3.6.2 👍 )

WARNING("CasADi was not compiled with WITH_OPENMP=ON. Falling back to serial evaluation.") [D:\bld\casadi_1682498776472\work\casadi\core\map.cpp:404]

[jgillis]

Nowadays, it's probably easiest to do your Python CasADi source builds through pip: https://github.com/casadi/casadi/wiki/FAQ:-how-do-I-configure-a-pip-source-build%3F

'thread' not working is news to me. Did you use our pip wheels, or the non-official conda package?

[halvorak]

First, I tried the conda-forge package.

Now, I also tried to make a new venv (in anaconda) and then install casadi separately by pip. That gave the same result as installing by conda-forge. Specifically, I ran this in anaconda prompt:

conda create -n caEnv python=3.10 scipy numpy pandas seaborn matplotlib scikit-learn statsmodels arviz scikit-image spyder salib jupyter
activate caEnv
pip install casadi

NB: I have changed my computer since my original question. I made a new venv with casadi 3.5.5, and now thread does not run in parallell there either.. So I guess the issue is somehow connected with my new pc (on my old PC with casadi 3.5.5, the %timeit sol = F_map_t(par_samples), where F_map_t is the threaded map, was approximately twice as fast as the other map-executions. On my new PC, it is actually slightly slower)

The sys.version and ca.version on my new computer gives the following. My old computer is broken, so cannot give sys.version there.
sys.version='3.9.16 (main, Mar 8 2023, 10:39:24) [MSC v.1916 64 bit (AMD64)]'
ca.version='3.5.5'

[jgillis]

'thread' seems to work fine on my windows machine.
Can you please provide the output of #3184 (comment) ?

[halvorak]

Thanks, the output from your script in #3184 is:
Time serial : 0.10583329200744629
Time thread : 0.02234482765197754
Speedup: 4.736366449355001 - expected : 8

So it seems like thread is working also for me!

I think the issue I experienced was that I added to many threads to my problem, as discussed in this Stackoverflow question: why-using-more-threads-makes-it-slower-than-using-less-threads. My new computer has 16 cpu-cores, but the fastest evaluation in the example I gave above is to use 4 threads. If I modify the code in my original question as below:

#define map functions
n_thread = 16
F_map = F2.map(N)
F_map_t = F2.map(N, "thread", n_thread)

#run time
print(f"F_map, {N=}, {n_thread=}, {os.cpu_count()=}")
%timeit sol = F_map(par_samples)
print(f"F_map_t")
%timeit sol = F_map_t(par_samples)

it prints out:
F_map, N=10000, n_thread=16, os.cpu_count()=16
429 ms ± 5.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
F_map_t
468 ms ± 6.46 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

If I change n_thread to 4, the I get:
F_map, N=10000, n_thread=4, os.cpu_count()=16
429 ms ± 3.28 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
F_map_t
305 ms ± 1.11 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

So thread was working for me after all, but I guess I did a rookie mistake by thinking more threads are always better. Anyhow, thanks a lot for the help, it was very useful!

PS: I think the script in #3184 is not completely correct, when I print out Y_ it is only the first row which is filled with values, the remaining rows are zeros.