Quasi-Random Sequence Example — faure_sequence, random_uniform, rng_get_seed

This example compares pseudo-random sampling via utilities.random_uniform() against a quasi-random (low-discrepancy) Faure/Halton sequence via utilities.faure_sequence() for 2-D point sets of 500 samples each. It also demonstrates utilities.rng_get_seed() to read back the active seed.

The approximate star discrepancy is used as a quantitative measure of space-filling quality: lower values indicate more uniform coverage of the unit square — a desirable property for numerical integration and design-of-experiments applications.

Example Code

"""IMSL quasi-random example: faure_sequence, random_uniform, rng_get_seed.

Outputs:
- Table printed to stdout
- SVG scatter plot saved to test_output/example_imsl_quasi_random.svg
"""
from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from utilities import faure_sequence, random_uniform, rng_get_seed, rng_set_seed


def _star_discrepancy_2d(pts: np.ndarray, n_grid: int = 20) -> float:
    """Approximate star discrepancy for a 2-D point set.

    Evaluates the maximum difference between the empirical proportion in
    [0, u] x [0, v] and the area u*v over a coarse grid.

    Args:
        pts (np.ndarray): Array of shape ``(n, 2)`` with points in ``[0,1]^2``.
        n_grid (int): Number of grid steps per dimension. Default is 20.

    Returns:
        float: Approximate star discrepancy (lower is better).
    """
    n = len(pts)
    max_disc = 0.0
    grid = np.linspace(0.0, 1.0, n_grid + 1)[1:]
    for u in grid:
        for v in grid:
            count = np.sum((pts[:, 0] <= u) & (pts[:, 1] <= v))
            disc = abs(count / n - u * v)
            if disc > max_disc:
                max_disc = disc
    return max_disc


def run_demo_imsl_quasi_random() -> Dict[str, object]:
    """Demonstrate quasi-random vs pseudo-random point generation.

    Generates 500 2-D points using both ``random_uniform`` (pseudo-random)
    and ``faure_sequence`` (quasi-random / low-discrepancy), then compares
    their space-filling properties via approximate star discrepancy.
    Also shows ``rng_get_seed`` for saving and restoring RNG state.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``pseudo_pts``,
            ``quasi_pts``, ``pseudo_discrepancy``, ``quasi_discrepancy``,
            ``seed_before``, ``seed_after``, and ``plot_path``.
    """
    n_pts = 500

    # ── rng_get_seed / rng_set_seed round-trip ────────────────────────────────
    rng_set_seed(2024)
    seed_before = rng_get_seed()
    print(f"\nRNG seed (before generation): {seed_before}")

    # ── Pseudo-random 2-D points ──────────────────────────────────────────────
    x_pseudo = random_uniform(n_pts, seed=2024)
    y_pseudo = random_uniform(n_pts, seed=9999)
    pseudo_pts = np.column_stack([x_pseudo, y_pseudo])

    seed_after = rng_get_seed()
    print(f"RNG seed (after random_uniform calls): {seed_after}")

    # ── Quasi-random Faure/Halton 2-D sequence ────────────────────────────────
    quasi_pts = faure_sequence(n_points=n_pts, n_dim=2, seed=None)

    # ── Discrepancy statistics ────────────────────────────────────────────────
    pseudo_disc = _star_discrepancy_2d(pseudo_pts)
    quasi_disc = _star_discrepancy_2d(quasi_pts)

    print("\nIMSL Quasi-Random Sequence Example")
    print("=" * 55)
    print(f"\n{'Method':<25} {'Points':>8} {'Star discrepancy':>18}")
    print("-" * 55)
    print(f"{'random_uniform (pseudo)':<25} {n_pts:>8} {pseudo_disc:>18.6f}")
    print(f"{'faure_sequence (quasi)':<25} {n_pts:>8} {quasi_disc:>18.6f}")
    print("-" * 55)
    print(
        f"\nQuasi-random discrepancy is "
        f"{pseudo_disc / quasi_disc:.1f}× lower than pseudo-random "
        f"(better space filling)."
    )

    # Sample first few points
    print("\nFirst 5 pseudo-random 2-D points:")
    for i in range(5):
        print(f"  [{pseudo_pts[i, 0]:.4f}, {pseudo_pts[i, 1]:.4f}]")

    print("\nFirst 5 quasi-random (Faure/Halton) 2-D points:")
    for i in range(5):
        print(f"  [{quasi_pts[i, 0]:.4f}, {quasi_pts[i, 1]:.4f}]")

    # ── Plot ──────────────────────────────────────────────────────────────────
    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_quasi_random.svg"

    fig, axes = plt.subplots(1, 2, figsize=(13, 6))

    ax1, ax2 = axes

    ax1.scatter(
        pseudo_pts[:, 0], pseudo_pts[:, 1],
        s=4, alpha=0.5, color="#b45309",
    )
    ax1.set_xlim(0, 1)
    ax1.set_ylim(0, 1)
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.set_title(
        f"random_uniform (pseudo-random)\n"
        f"n={n_pts}  discrepancy={pseudo_disc:.4f}"
    )
    ax1.grid(True, alpha=0.2)

    ax2.scatter(
        quasi_pts[:, 0], quasi_pts[:, 1],
        s=4, alpha=0.6, color="#1d4ed8",
    )
    ax2.set_xlim(0, 1)
    ax2.set_ylim(0, 1)
    ax2.set_xlabel("x")
    ax2.set_ylabel("y")
    ax2.set_title(
        f"faure_sequence (quasi-random)\n"
        f"n={n_pts}  discrepancy={quasi_disc:.4f}"
    )
    ax2.grid(True, alpha=0.2)

    fig.suptitle(
        "IMSL Quasi-Random: random_uniform vs faure_sequence",
        fontweight="bold",
    )
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    print(f"\nPlot saved to: {plot_path}")

    return {
        "pseudo_pts": pseudo_pts,
        "quasi_pts": quasi_pts,
        "pseudo_discrepancy": pseudo_disc,
        "quasi_discrepancy": quasi_disc,
        "seed_before": seed_before,
        "seed_after": seed_after,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_quasi_random()

Plot Output

Side-by-side scatter plots of pseudo-random vs quasi-random 2-D points in the unit square

Left: 500 pseudo-random points from random_uniform showing visible clustering. Right: 500 Faure/Halton quasi-random points with visibly more uniform coverage. Star discrepancy values are annotated in each panel title.

Console Output

RNG seed (before generation): 2024
RNG seed (after random_uniform calls): 9999

IMSL Quasi-Random Sequence Example
=======================================================

Method                      Points   Star discrepancy
-------------------------------------------------------
random_uniform (pseudo)        500           0.056500
faure_sequence (quasi)         500           0.010000
-------------------------------------------------------

Quasi-random discrepancy is 5.7× lower than pseudo-random (better space filling).

First 5 pseudo-random 2-D points:
  [0.6758, 0.8387]
  [0.2143, 0.7620]
  [0.3095, 0.2624]
  [0.7995, 0.7889]
  [0.9958, 0.2437]

First 5 quasi-random (Faure/Halton) 2-D points:
  [0.0000, 0.0000]
  [0.5000, 0.3333]
  [0.2500, 0.6667]
  [0.7500, 0.1111]
  [0.1250, 0.4444]

Plot saved to: test_output\example_imsl_quasi_random.svg