<p dir="auto">From <a href="https://stackoverflow.com/questions/68927829/how-to-use-use-numpy-random-choice-to-get-progressively-longer-sequences-with-th" rel="nofollow">https://stackoverflow.com/questions/68927829/how-to-use-use-numpy-random-choice-to-get-progressively-longer-sequences-with-th</a></p>
<p dir="auto">When generating random numbers with <code class="notranslate">np.random.Generator.choice</code> and having probabilities, the chosen elements in the first positions will depend on the number of elements to generate, leading to inconsistent results if one tries to replicate results but with varying input lengths:</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="import numpy as np
n = 100
prob = np.random.choice(100, n)
prob = prob / np.sum(prob)
for i in np.arange(10, 40, 10):
    rng = np.random.default_rng(seed=123)     
    selected = rng.choice(n, size=i, replace=False, p=prob)
    print(str(i) + ': ' + str(selected))"><pre class="notranslate"><span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>
<span class="pl-s1">n</span> <span class="pl-c1">=</span> <span class="pl-c1">100</span>
<span class="pl-s1">prob</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-s1">random</span>.<span class="pl-en">choice</span>(<span class="pl-c1">100</span>, <span class="pl-s1">n</span>)
<span class="pl-s1">prob</span> <span class="pl-c1">=</span> <span class="pl-s1">prob</span> <span class="pl-c1">/</span> <span class="pl-s1">np</span>.<span class="pl-en">sum</span>(<span class="pl-s1">prob</span>)
<span class="pl-k">for</span> <span class="pl-s1">i</span> <span class="pl-c1">in</span> <span class="pl-s1">np</span>.<span class="pl-en">arange</span>(<span class="pl-c1">10</span>, <span class="pl-c1">40</span>, <span class="pl-c1">10</span>):
    <span class="pl-s1">rng</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-s1">random</span>.<span class="pl-en">default_rng</span>(<span class="pl-s1">seed</span><span class="pl-c1">=</span><span class="pl-c1">123</span>)     
    <span class="pl-s1">selected</span> <span class="pl-c1">=</span> <span class="pl-s1">rng</span>.<span class="pl-en">choice</span>(<span class="pl-s1">n</span>, <span class="pl-s1">size</span><span class="pl-c1">=</span><span class="pl-s1">i</span>, <span class="pl-s1">replace</span><span class="pl-c1">=</span><span class="pl-c1">False</span>, <span class="pl-s1">p</span><span class="pl-c1">=</span><span class="pl-s1">prob</span>)
    <span class="pl-en">print</span>(<span class="pl-en">str</span>(<span class="pl-s1">i</span>) <span class="pl-c1">+</span> <span class="pl-s">': '</span> <span class="pl-c1">+</span> <span class="pl-en">str</span>(<span class="pl-s1">selected</span>))</pre></div>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="10: [71  5 20 16 81 92 26 82 88 52]
20: [71  5 20 16 81 92 26 82 88 52 24 19 75 65 22 80 53 25 13 48]
30: [71  5 20 16 81 92 26 82 88 52 24 19 75 65 22 80 53 15 50 60  1 48 74 91 94 87 18 73 31 79]"><pre class="notranslate"><code class="notranslate">10: [71  5 20 16 81 92 26 82 88 52]
20: [71  5 20 16 81 92 26 82 88 52 24 19 75 65 22 80 53 25 13 48]
30: [71  5 20 16 81 92 26 82 88 52 24 19 75 65 22 80 53 15 50 60  1 48 74 91 94 87 18 73 31 79]
</code></pre></div>
<p dir="auto">Which is quite unexpected and does not happen in e.g. R's <code class="notranslate">sample</code>.</p>
<p dir="auto">I would assume the optimal algorithm for drawing these random numbers amounts to having some binary tree saving weights of its left and right subtrees, and then drawing random numbers in the tree's range and passing down the tree until reaching the bottom and removing that weight, which should not make any choice at position <code class="notranslate">n</code> depend on what's going to happen at position <code class="notranslate">n+1</code>.</p>
<p dir="auto">In the FFTPACK code, we have some "clever" code that uses the pattern <code class="notranslate">NPY_SIGINT_ON</code>, then run some inner loop, then <code class="notranslate">NPY_SIGINT_OFF</code>. And what this does is:</p>
<p dir="auto"><code class="notranslate">NPY_SIGINT_ON</code> uses <code class="notranslate">setjmp</code> to make a jump buffer target, saves the current signal handler for SIGINT (control-C), and then registers a new SIGINT handler that attempts to handle control-C by <code class="notranslate">longjmp</code>ing out to the <code class="notranslate">setjmp</code> buffer. The idea is that this allows us to cancel long-running calculations.</p>
<p dir="auto">The actual implementations of these macros is:</p>
<div class="highlight highlight-source-c notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="#    define NPY_SIGINT_ON {                                             \
                   PyOS_sighandler_t _npy_sig_save;                     \
                   _npy_sig_save = PyOS_setsig(SIGINT, _PyArray_SigintHandler);\
 \
                   if (NPY_SIGSETJMP(*((NPY_SIGJMP_BUF *)_PyArray_GetSigintBuf(\
)), \
                                 1) == 0) {

#    define NPY_SIGINT_OFF }                                      \
        PyOS_setsig(SIGINT, _npy_sig_save);                       \
        }"><pre class="notranslate"><span class="pl-k">#    define</span> <span class="pl-c1">NPY_SIGINT_ON</span> {                                             \
                   PyOS_sighandler_t _npy_sig_save;                     \
                   _npy_sig_save = PyOS_setsig(SIGINT, _PyArray_SigintHandler);\
 \
                   if (NPY_SIGSETJMP(*((NPY_SIGJMP_BUF *)_PyArray_GetSigintBuf(\
)), \
                                 1) == 0) {

<span class="pl-k">#    define</span> <span class="pl-c1">NPY_SIGINT_OFF</span> }                                      \
        PyOS_setsig(SIGINT, _npy_sig_save);                       \
        }</pre></div>
<p dir="auto">So this has two race conditions:</p>
<p dir="auto">Minor problem: If we receive a control-C <em>in between</em> the call to <code class="notranslate">PyOS_setsig</code> and the call to <code class="notranslate">NPY_SIGSETJMP</code>, then we'll longjmp out of the signal handler into a uninitialized buffer, and <a href="https://xkcd.com/1133/" rel="nofollow">will not go to space today</a>.</p>
<p dir="auto">Major problem: suppose that the following sequence occurs:</p>
<ol dir="auto">
<li>thread 1 enters <code class="notranslate">NPY_SIGINT_ON</code>, stashes Python's default sigint handler in its <code class="notranslate">_npy_sig_save</code></li>
<li>thread 2 enters <code class="notranslate">NPY_SIGINT_ON</code>, stashes <em>our</em> sigint handler in its <code class="notranslate">_npy_sig_save</code></li>
<li>thread 1 leaves via <code class="notranslate">NPY_SIGINT_OFF</code>, restores Python's default sigint handler from its <code class="notranslate">_npy_sig_save</code></li>
<li>thread 2 leaves via <code class="notranslate">NPY_SIGINT_OFF</code>, then restores <em>our</em> sigint handler from its <code class="notranslate">_npy_sig_save</code></li>
</ol>
<p dir="auto">Now our signal handler gets left installed indefinitely, and eventually when control-C gets hit it will attempt to jump out to a stack frame that was reused ages ago. This is a bad problem and you will not go to space today.</p>
<p dir="auto">By adding some <code class="notranslate">printf</code> checks to <code class="notranslate">NPY_SIGINT_OFF</code>, I've confirmed that this is the cause of the segfault that <a href="https://mail.scipy.org/pipermail/numpy-discussion/2016-April/075361.html" rel="nofollow">Oscar reported on numpy-discussion</a>. It's very easy to replicate even from Python code:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import numpy
numpy.test()
import os
import signal
os.kill(os.getpid(), signal.SIGINT)
# -&gt; segfault"><pre class="notranslate"><code class="notranslate">import numpy
numpy.test()
import os
import signal
os.kill(os.getpid(), signal.SIGINT)
# -&gt; segfault
</code></pre></div>
<p dir="auto">Also, the whole architecture here is somewhat screwed up, because the Python model is that all signals are directed to the main thread (by setting them to masked on all other threads). So it doesn't even make sense to be installing a signal handler to interrupt operations on a thread that can't receive signals.</p>
<p dir="auto">(If you look at the definition of <code class="notranslate">_PyArray_SigintHandler</code> and <code class="notranslate">_PyArray_GetSigintBuf</code> in <code class="notranslate">multiarraymodule.c</code>, then you'll see that there's some basic defenses against this: we keep our jump buffer in TLS, and also a flag <code class="notranslate">sigint_buf_init</code>, so that idea is that if we start running the signal handler in a thread that hasn't actually allocated a jump buffer, then we throw away the signal. This is also pretty buggy: (1) nothing ever sets the flag back to 0, so if a thread has <em>ever</em> had a valid jump buffer on its stack then we will happily jump to it, even if it's no longer valid. (2) we shouldn't throw away the signal, we should pass it on to Python.)</p>
<p dir="auto">So the simple solution is (1) check whether we are actually the one thread that can receive signals, (2) if so, then register our handler etc.; (3) if not, then make these macros into no-ops.</p>
<p dir="auto">The one wrinkle here is that in theory it's possible for someone to call <code class="notranslate">sigprocmask</code> themselves and enable signal receiving on a non-default thread. We might want to be robust to this. A nice extra constraint is that there's no C API to check whether we are in the main thread; the only way I can think of to check that is actually to call <code class="notranslate">sigprocmask</code> and check whether SIGINT delivery is enabled for this thread.</p>