<div dir="ltr"><div style="font-family:verdana,sans-serif;font-size:small" class="gmail_default"></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">Hi all,</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">I have been working on integrating UNU.RAN in scipy.stats for a few weeks on my fork (<a href="https://github.com/tirthasheshpatel/scipy/pull/6" target="_blank">https://github.com/tirthasheshpatel/scipy/pull/6</a>) with my mentors, Christoph and Nicholas. The design is still fairly basic but is almost ready for initial reviews. Also, it builds on Windows, macOS, and Linux. Tests pass with an exception of unrelated failures on Linux. I am planning to propose a PR on SciPy this week but before I do that, I thought it would be a good idea to summarize the progress here.</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">The methods I propose to add and details regarding the parameters they take can be found in this excel sheet: <a href="https://docs.google.com/spreadsheets/d/1r36HypXwpit7sHt9YAe3K7yPL7HNNuQjTkbVYQIJ4xI/edit?usp=sharing" target="_blank">https://docs.google.com/spreadsheets/d/1r36HypXwpit7sHt9YAe3K7yPL7HNNuQjTkbVYQIJ4xI/edit?usp=sharing</a>. Feel free to comment on the sheet if you have any suggestions.</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">The API looks something like this:</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    import numpy as np<br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    from scipy.stats import TransformedDensityRejection</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    pdf = lambda x: 1 - x*x</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    dpdf = lambda x: -2*x</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    urng = np.random.default_rng(123)</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    rng = TransformedDensityRejection(pdf, dpdf, domain=(-1,1), seed=urng)</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    # sample a 100,000 samples from the distribution<br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    rvs = rng.rvs(100_000)<br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    # evaluate the upper bound of the PPF at some points in the domain</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    u = np.linspace(0, 1, num=1000)</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">    ppf_hat = rng.ppf_hat(u)</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">Internal API</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">------------<br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">For Python callbacks, I have used the `ccallback` API and used setjmp/longjmp for error handling. Using `ctypes` gives similar performance but I wasn't able to figure out a way to handle errors occurring in UNU.RAN and in the callbacks. Here is the C code I wrote to acquire and release the callbacks: <a href="https://github.com/tirthasheshpatel/scipy/blob/unuran_c_thunk/scipy/stats/unuran/unuran_callback.h" target="_blank">https://github.com/tirthasheshpatel/scipy/blob/unuran_c_thunk/scipy/stats/unuran/unuran_callback.h</a></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">For more details on the internal design of the API, I have attempted to create a few flowcharts on <a href="http://draw.io" target="_blank">draw.io</a> but they are very much a WIP: <a href="https://drive.google.com/file/d/1TH70SSvvc5eF6-YmDO8kFNvLqKaGROW-/view?usp=sharing" target="_blank">https://drive.google.com/file/d/1TH70SSvvc5eF6-YmDO8kFNvLqKaGROW-/view?usp=sharing</a></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">Other Notes</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">-----------<br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">I think NumPy 1.19 introduced a stable interface for its Generator API. Before that the RandomState API was stable. But RandomState doesn't have a Cython Interface and its Python interface is very slow to be practically useful. For instance, it takes around 7ms to sample 100,000 samples from the normal distribution using the TDR method on NumPy >= 1.19 while it takes around 911ms on NumPy == 1.18 which is around 130 times slower! Is there a plan to drop 1.16 soon and can we use the unstable Generator API from 1.17 onwards or would it too unsafe? Maybe this discussion isn't suited here but I just thought to put it out as a note.</div></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">Anyways, please feel free to comment on the design, naming convention, etc and suggest changes/ideas that would be good to incorporate before proposing a PR on SciPy. Thanks for reading this!<br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">--<br clear="all"></div><div><div dir="ltr" data-smartmail="gmail_signature"><div dir="ltr"><div>Kind Regards,</div><div>Tirth Patel</div></div></div></div></div>