WCS conversion functions modify input array in-place, making them thread-unsafe

[astrofrog]

While investigating multi-threaded issues in reproject, I came across an issue which is that if the input array to e.g. wcs_pix2world is a (N, 2) array and if the origin is set to 0, the input array is modified in-place temporarily to make the pixel values be 1-based since this is what WCSLIB expects. This causes multi-threading to fail spectacularly in this case if multiple threads use the same input array, as demonstrated in the following example:

import numpy as np
from astropy.wcs import WCS
from multiprocessing.pool import ThreadPool

wcs = WCS(naxis=2)

N = 1_000_000

pixel = np.random.randint(-1000, 1000, N * 2).reshape((N, 2)).astype(float)


def round_trip_transform(pixel):
    world = wcs.wcs_pix2world(pixel, 0)
    pixel = wcs.wcs_world2pix(world, 0)
    return pixel


for n_threads in [1, 2, 8]:

    print("---")
    print("N_threads:", n_threads)

    pool = ThreadPool(n_threads)
    results = pool.map(round_trip_transform, (pixel,) * n_threads)

    for pixel2 in results:
        print(f"{np.sum(~np.isclose(pixel, pixel2))} {np.unique(pixel - pixel2)}")

The output in this case is:

---
N_threads: 1
0 [0.]
---
N_threads: 2
2000000 [-1.]
2000000 [-1.]
---
N_threads: 8
2000000 [-5. -4. -3.]
2000000 [-5. -4.]
2000000 [-5. -4.]
2000000 [-7. -6. -5.]
2000000 [-7. -6. -5.]
2000000 [-7. -6. -5.]
2000000 [-5. -4. -3.]
2000000 [-5. -4.]

What's happening here is that several threads are applying an offset in-place in the input array at the same time:

astropy/astropy/wcs/src/wcslib_wrap.c

Line 1413 in 01e8729

preoffset_array(pixcrd, origin);

Note that if the origin is changed to 1, there is no issue since no offset is applied:

    world = wcs.wcs_pix2world(pixel, 1)
    pixel = wcs.wcs_world2pix(world, 1)

gives:

---
N_threads: 1
0 [0.]
---
N_threads: 2
0 [0.]
0 [0.]
---
N_threads: 8
0 [0.]
0 [0.]
0 [0.]
0 [0.]
0 [0.]
0 [0.]
0 [0.]
0 [0.]

An easy fix is to copy the array before passing it to the C extension, though obviously this will increase memory usage, which isn't ideal. But in any case we shouldn't be modifying the input array in-place. We should probably check if WCSLIB now has a mechanism for dealing with different origins on-the-fly during the conversions and whether our code is still needed.

[pllim]

So, is this related to #16245 ?

[astrofrog]

No, these are two separate issues related to astropy.wcs and multi-threading, with different causes.

[manodeep]

Noting here that my proposed pthreads solution for #16245 to ensure a single-threaded execution likely will fix this issue. But that solution does require adding a compile and runtime dependency for pthreads.

Essentially consists of using pthread_mutex_trylock to only execute the single-threaded region on the thread that gets the lock, and then pthread_cond_wait for the other threads to wait while the single thread finishes executing (i.e., the calls to pre-offset and un-offset)