Make the three examples runnable on the modernized stack

End-to-end verification of `examples/` revealed three issues, none of
them regressions in the library itself but all of them rough edges
that left the examples broken on a clean v1.0.0 install:

1. `examples/wlsqm_example.py` and `examples/lapackdrivers_example.py`
   both used `plt.grid(b=True, which='both')`, where `b=` is the old
   matplotlib 2.x kwarg. Matplotlib 3.x renamed it to `visible=` and
   now errors out with a confusing "keyword grid_b is not recognized".
   Changed to `plt.grid(visible=True, which='both')` at all three sites.

2. `examples/wlsqm_example.py` imports sympy to symbolically
   differentiate a manufactured solution and check the wlsqm fit
   against the analytical derivatives. sympy was not in the dev deps,
   so the example failed at import time. Added sympy to the dev list
   in pyproject.toml. Also restored the comments PDM stripped from
   the dev list during the previous `pdm add --dev` invocation
   (build, pyyaml, meson-python all had explanatory comments that
   went missing).

3. `examples/lapackdrivers_example.py` had a pre-existing test-logic
   flaw: it asserted `|x_wlsqm - x_numpy| < 1e-10` against a
   `numpy.linalg.solve` reference, treating NumPy as ground truth.
   On this particular run, wlsqm and `scipy.linalg.solve` agree
   bit-identically (they call the same LAPACK DGESV/DSYSV) but both
   differ from NumPy by ~1.5e-6. Investigation: residuals
   `‖A x − b‖ ≈ 1e-13` for all three. The matrices are `(U + U.T)/2`
   for U ~ uniform[0,1], moderately ill-conditioned (κ ~ 1e4 at n=117),
   indefinite, and the three LAPACK paths legitimately pick slightly
   different equally-valid solutions — that is a property of finite-
   precision arithmetic on a borderline-conditioned matrix, not a bug.

   The example never tripped consistently in 2017 because it used
   unseeded `np.random.sample(...)` and depended on lucky draws.

   Replaced the per-element `|x_wlsqm - x_numpy|` assertion with a
   per-solver relative residual check: `‖A x − b‖ / ‖b‖ < 1e-8`,
   computed per problem instance with a vectorized
   `np.einsum('ijk,jk->ik', A, x) - b`. This is the right sanity check
   for a linear solver — it does not depend on having a "ground truth"
   solution. The 1e-8 threshold accommodates DSYSV (Bunch-Kaufman)
   producing noticeably larger residuals than DGESV on indefinite
   matrices, which is normal and expected. It is still 8 orders of
   magnitude above machine epsilon and tight enough to catch any
   realistic regression in the Cython wrappers.

   Also seeded `np.random.seed(42)` in both examples' `main()` so the
   per-run residuals (and the printed "max error" stats in
   wlsqm_example.py) are reproducible across runs and machines.
   Without seeding, the lapackdrivers residuals swung from ~1e-11 to
   ~1e-9 between runs purely as a function of which random matrix the
   global RNG drew.

After this commit, all three examples run cleanly to exit code 0 with
matplotlib in headless mode (`MPLBACKEND=Agg`):

  - examples/expertsolver_example.py: ~1 s
  - examples/lapackdrivers_example.py: ~30 s (lots of timed solver runs)
  - examples/wlsqm_example.py: ~1 min (test3d, test2d, test1d, testmany2d)

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: Technologicat

examples/lapackdrivers_example.py

@@ -65,6 +65,13 @@ def idfun(x): return x
 
 
 def main():
+    # Seed the legacy global RNG used by `np.random.sample(...)` calls below
+    # so that the residual-check thresholds are reproducible across runs and
+    # across machines. Without this, each run draws different random
+    # matrices, the conditioning varies wildly, and the per-solver residuals
+    # can swing by orders of magnitude run-to-run.
+    np.random.seed(42)
+
 #    # exact solution is (3/10, 2/5, 0)
 #    A = np.array( ( (2., 1.,  3.),
 #                    (2., 6.,  8.),
@@ -246,15 +253,38 @@ def main():
         drivers.mgeneralfactoredp( fact, ipiv, x5, ntasks )
         results7[j] = (time.time() - t0) / r
 
-        if use_numpy:
-#            print( np.max(np.abs(x - x3)) )  # DEBUG
-#            print( np.max(np.abs(x - x5)) )  # DEBUG
-            print( np.max(np.abs(x2 - x4)) )  # DEBUG
-            assert (np.abs(x - x5) < 1e-10).all(), "Something went wrong, solutions do not match"  # check general solver first
-            assert (np.abs(x - x3) < 1e-10).all(), "Something went wrong, solutions do not match"  # check general solver
-#            assert (np.abs(x - x2) < 1e-5).all(), "Something went wrong, solutions do not match"  # doesn't make sense to compare, DSYSV is more accurate for badly conditioned symmetric matrices
-            assert (np.abs(x2 - x4) < 1e-7).all(), "Something went wrong, solutions do not match"  # check symmetric solvers against each other
-                                                                                                   # (not exactly the same algorithm (DSYTRS2 vs. DSYTRS), so there may be slight deviation)
+        # Verify each solver produced a valid solution.
+        #
+        # The right sanity check for a linear solver is the relative residual
+        # ‖A x − b‖ / ‖b‖, NOT ‖x − x_reference‖. Two valid LAPACK calls on a
+        # moderately ill-conditioned matrix can produce solutions that differ
+        # at, say, 1e-6 — that is a property of the conditioning, not a bug —
+        # while both still have a tiny residual ‖A x − b‖ ~ machine epsilon.
+        # This matters here because msymmetrizep produces matrices with
+        # κ(A) ~ 1e4 at n=117, and the historical (vs-NumPy) check used to
+        # trip nondeterministically as a function of the unseeded RNG.
+        #
+        # Threshold rationale: 1e-8 covers both DGESV and DSYSV across the
+        # whole size range used here (max n ~ 117). DSYSV (Bunch-Kaufman)
+        # tends to produce noticeably larger residuals than DGESV on
+        # indefinite matrices because of its different pivoting strategy,
+        # and the random `(U + U.T) / 2` matrices we generate here are
+        # almost always indefinite and moderately ill-conditioned. 1e-8 is
+        # still 8 orders of magnitude above machine epsilon and tight
+        # enough to catch any realistic regression in the Cython wrappers.
+        #
+        # einsum 'ijk,jk->ik' computes A[:,:,k] @ x[:,k] for each problem
+        # instance k, vectorized.
+        b_norm = np.linalg.norm(b, axis=0)
+        b_norm = np.maximum(b_norm, 1.0)  # guard against the all-zero RHS edge case
+        for label, solver_x in (("msymmetricp",                x2),
+                                ("mgeneralp",                  x3),
+                                ("msymmetricfactorp+factored", x4),
+                                ("mgeneralfactorp+factored",   x5)):
+            residual = np.linalg.norm(np.einsum('ijk,jk->ik', A, solver_x) - b, axis=0) / b_norm
+            worst = residual.max()
+            assert worst < 1e-8, f"{label}: relative residual {worst:.3e} exceeds 1e-8"
+            print(f"        {label} max relative residual: {worst:.3e}")
 
 
 # old, serial only
@@ -304,7 +334,7 @@ def main():
     plt.ylabel('t')
     plt.title('Average time per problem instance, %d parallel tasks' % (ntasks))
     plt.axis('tight')
-    plt.grid(b=True, which='both')
+    plt.grid(visible=True, which='both')
     plt.legend(loc='best')
 
     plt.savefig('figure1_latest.pdf')

examples/wlsqm_example.py

@@ -899,7 +899,7 @@ def test2d():
     ax.plot( (xi[0],), (xi[1],), linestyle='none', marker='x', markeredgecolor='k', markerfacecolor='none' )
     plt.axis('tight')
     axis_marginize(ax, 0.02, 0.02)
-    plt.grid(b=True, which='both')
+    plt.grid(visible=True, which='both')
     plt.xlabel('x')
     plt.ylabel('y')
     plt.subplot(1,2, 2)
@@ -1256,12 +1256,20 @@ def test1d():
 
     plt.axis('tight')
     axis_marginize(ax, 0.02, 0.02)
-    plt.grid(b=True, which='both')
+    plt.grid(visible=True, which='both')
     plt.xlabel('x')
     plt.ylabel('y')
 
 
 def main():
+    # Seed the legacy global RNG used by `np.random.sample(...)` and
+    # `np.random.normal(...)` calls in the test functions below, so that
+    # each run produces the same point clouds, the same noise realizations,
+    # and the same printed error figures. This makes the example
+    # reproducible across runs and across machines, which matters because
+    # the printed "max error" lines are how a reader gauges fit quality.
+    np.random.seed(42)
+
     test3d()
     test2d()
     test1d()

pyproject.toml

@@ -157,6 +157,10 @@ dev = [
     "jedi>=0.19.2",
     "scipy>=1.9",
     "matplotlib",
+    # examples/wlsqm_example.py uses sympy to symbolically differentiate a
+    # manufactured solution and check the wlsqm fit against the analytical
+    # derivatives. Not used by the library itself or by the test suite.
+    "sympy>=1.14.0",
     # Needed in the venv (not just the isolated build env) so meson-python's
     # editable loader can rebuild the extension on import after .pyx/.pxd edits.
     "meson-python>=0.17",