Sobol & Halton Generators

[tensor-calculus]

Implements part of #209

Added Sobol and Halton quasi-random low-discrepancy sequence generators.

I am not very familiar with QMC methods, so the implementations were largely based on the structure of the existing generators in random.py and my understanding of the SciPy documentation. I also created test_random.py and added a few basic tests covering output shape and reproducibility.

I am using rejection sampling for this which could be a potential concern and might need additional tests. If my understanding is correct, cases where the hull occupies only a small fraction of its bounding box (for example, a thin diagonally oriented polygon) could result in poor acceptance rates which would require many iterations before the requested number of points is generated and the current implementation does not impose a maximum iteration limit, which means it could go on for a long time.

If the overall approach looks reasonable, I can continue working on the remaining generators. From what I could gather, their implementation will also be largely similar.

[codecov]

Codecov Report

❌ Patch coverage is 64.17910% with 24 lines in your changes missing coverage. Please review.
✅ Project coverage is 73.2%. Comparing base (a798e9e) to head (5a3a7b2).
⚠️ Report is 5 commits behind head on main.

Files with missing lines	Patch %	Lines
pointpats/random.py	64.2%	24 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##            main    #210     +/-   ##
=======================================
+ Coverage   72.9%   73.2%   +0.3%     
=======================================
  Files         12      12             
  Lines       2118    2178     +60     
=======================================
+ Hits        1543    1594     +51     
- Misses       575     584      +9     

Files with missing lines	Coverage Δ
pointpats/random.py	67.0% <64.2%> (-0.9%)	⬇️

... and 2 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ljwolf]

The implementation looks great---simple, clean and in style. Thanks! Could you make a simple ipython notebook for them?

You're right that the rejection sampling will be inefficient for certain geometries (arch/banana-shaped ones as well as thin diagonal ones). We can solve this somewhat by sampling chunks of points in one shot, but this can be hard to tune correctly. I went with single-unit sampling in other classes due to some tuning experiments that showed it's hard to come up with a useful default (i.e. sqrt(n) was often too big, log(n) too small, area(bbox)/area(shape)*n was effective for the edge cases but too large for the normal cases).

If you'd like, keep going for the R2? My implementation is over here, but has a different API/use case.

[sjsrey]

Should this (and the while block) be nested inside the for loop on 882?

[tensor-calculus]

OMG, Yes! My mistake. I copied this from the Sobol function I wrote earlier and must have accidentally messed up the indentation. I'll fix it ASAP.

I was wondering why the tests didn't catch it, but then I realized I was only asserting the shape of the result, which is hard-coded in the function. I'll need to improve the tests to validate the actual output as well.

[sjsrey]

This looks really good.
Only a few questions/issues.

[sjsrey]

I'm wondering if this can be vectorized?

[tensor-calculus]

I think it's possible, but it would require some changes to the contains() method. Currently it does:

@contains.register
def _(shape: numpy.ndarray, x: float, y: float):
    """
    If provided an ndarray, assume it's a bbox
    and return whether the point falls inside
    """
    xmin, xmax = shape[0], shape[2]
    ymin, ymax = shape[1], shape[3]
    in_x = (xmin <= x) and (x <= xmax)
    in_y = (ymin <= y) and (y <= ymax)
    return in_x & in_y

The implementation uses and operator so it won't work with array inputs. We'll have to switch to & operator to support vectorization. Until then, I think we'll have to evaluate points individually in a loop.

[tensor-calculus]

The implementation looks great---simple, clean and in style. Thanks! Could you make a simple ipython notebook for them?

You're right that the rejection sampling will be inefficient for certain geometries (arch/banana-shaped ones as well as thin diagonal ones). We can solve this somewhat by sampling chunks of points in one shot, but this can be hard to tune correctly. I went with single-unit sampling in other classes due to some tuning experiments that showed it's hard to come up with a useful default (i.e. sqrt(n) was often too big, log(n) too small, area(bbox)/area(shape)*n was effective for the edge cases but too large for the normal cases).

If you'd like, keep going for the R2? My implementation is over here, but has a different API/use case.

Thanks! I'll look into the R2 implementation as well. Should I add the examples for these to the existing random.ipynb notebook?

[tensor-calculus]

Thanks! I'll look into the R2 implementation as well. Should I add the examples for these to the existing random.ipynb notebook?

I've added the R2 generator as well. It's largely similar to the previous generators, but please take a look and let me know if the logic looks correct. I'll add examples for all three shortly. Thanks!