Make Python optional: plotting and pyamg become package extensions; bump to 1.0.2

PyPlot and PyCall move from [deps] to [weakdeps]. The solver core now has no
Python dependency at all (the dependency closure shrinks from 206 to 174
packages); plotting and the pyamg prolongator are opt-in extensions.

MultiGridBarrierPyPlotExt (load PyPlot; PyCall arrives with it) carries every
plot method -- 2D FEM, 1D, spectral, MGBSOL, ParabolicSOL with the matplotlib
animation, and plot3d.jl wholesale with its lazy pyvista import. The methods
extend PyPlot.plot and PyPlot.savefig, so
gives exactly the old behavior; the parent no longer exports plot/savefig.
MGB3DFigure and HTML5anim hold plain PNG/HTML data and stay in the parent.

MultiGridBarrierPyAMGExt (load PyCall) supplies amg_pyamg; the parent keeps
the exported stub. FFMPEG/PNGFiles/Base64 stay hard deps (not Python; the
extension uses them through the parent).

Docs gain Plotting and PyAMG tabs; index/zoo examples show the PyPlot load;
README quickstart updated and gains a Python-optional bullet. CI configures
PyCall against the system Python in a throwaway environment (the build config
is depot-wide) and keeps the OpenSSL_jll 3.0 pin CI-only, in the package env,
where it always lived; the test target gains PyPlot/PyCall extras and
runtests.jl loads PyPlot.

Verified: full Pkg.test green (all five extensions, including pyamg through
the new extension against the real Python pyamg); local docs build green;
Python-free load and solve confirmed in a PyCall-free environment.

Author: Sébastien Loisel

.github/workflows/CI.yml

@@ -48,11 +48,20 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v3
-      - name: Build package with system Python
+      # PyCall/PyPlot are weak dependencies (test-only extras), so PyCall is
+      # configured against the system Python in a throwaway environment (the
+      # build config is depot-wide). The OpenSSL_jll pin goes into the package
+      # env so the test env inherits it: it keeps Julia's libssl compatible
+      # with the system Python's at runtime.
+      - name: Configure PyCall against system Python
         env:
           PYTHON: python
-        run: julia --project=. -e 'using Pkg; Pkg.add(name="OpenSSL_jll", version="3.0"); Pkg.instantiate(); Pkg.build("PyCall")'
+        run: julia -e 'using Pkg; Pkg.activate(temp=true); Pkg.add(Pkg.PackageSpec(name="OpenSSL_jll", version="3.0")); Pkg.add("PyCall"); Pkg.build("PyCall")'
+      - name: Instantiate
+        run: julia --project=. -e 'using Pkg; Pkg.add(name="OpenSSL_jll", version="3.0"); Pkg.instantiate()'
       - name: Run tests with XVFB
+        env:
+          PYTHON: python
         run: xvfb-run --auto-servernum julia --project=. -e 'using Pkg; Pkg.test(coverage=true)'
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v7

Project.toml

@@ -1,7 +1,7 @@
 name = "MultiGridBarrier"
 uuid = "9e2c1f1d-9131-4ad4-b32f-bd2a0b0ecd1e"
 authors = ["Sébastien Loisel"]
-version = "1.0.1"
+version = "1.0.2"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -11,8 +11,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PNGFiles = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
-PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
-PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 QuadratureRules = "a08977f5-a20d-4b99-b11f-f5ea535e0575"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -23,11 +21,15 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CUDSS_jll = "4889d778-9329-5762-9fec-0578a5d30366"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 
 [extensions]
 MultiGridBarrierCUDAExt = ["CUDA", "CUDSS_jll"]
 MultiGridBarrierGmshExt = "Gmsh"
 MultiGridBarrierJuMPExt = "JuMP"
+MultiGridBarrierPyAMGExt = "PyCall"
+MultiGridBarrierPyPlotExt = ["PyPlot", "PyCall"]
 
 [compat]
 AlgebraicMultigrid = "1"
@@ -51,7 +53,10 @@ CUDSS_jll = "4889d778-9329-5762-9fec-0578a5d30366"
 Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Coverage", "CUDA", "CUDSS_jll", "JuMP", "Gmsh"]
+# PyPlot/PyCall activate the plotting and pyamg extensions under test.
+test = ["Test", "Coverage", "CUDA", "CUDSS_jll", "JuMP", "Gmsh", "PyPlot", "PyCall"]

README.md

@@ -20,9 +20,9 @@ Chebyshev spectral discretizations, with optional CUDA GPU acceleration.
 ## Install and quickstart
 
 ```julia
-using Pkg; Pkg.add("MultiGridBarrier")
+using Pkg; Pkg.add(["MultiGridBarrier", "PyPlot"])
 
-using MultiGridBarrier
+using MultiGridBarrier, PyPlot                    # PyPlot enables the plotting extension
 geom = fem2d_P2()                                 # a 2D P2 triangular mesh
 sol  = mgb_solve(assemble(amg(geom); p = 1.0))    # solve a (nonsmooth) p = 1 problem
 plot(sol)
@@ -38,6 +38,9 @@ plot(sol)
 - **Discretizations:** finite elements in 1D/2D/3D (simplicial `P1`/`P2` and tensor-product `Q_k`),
   plus Chebyshev spectral elements; all isoparametric.
 - **Solver:** an algebraic-multigrid hierarchy (`amg`) driving a barrier (interior-point) method.
+- **Python-optional:** the solver core has no Python dependency. Plotting (`using PyPlot`:
+  matplotlib in 1D/2D, PyVista in 3D) and the `pyamg` prolongators (`using PyCall`) are
+  opt-in extensions.
 - **JuMP modeling front end** (loads with `using JuMP`): state variational problems with
   `@variable`/`@constraint`/`@objective`; `optimize!` lowers them straight to the multigrid
   barrier pipeline.

docs/make.jl

@@ -25,8 +25,10 @@ makedocs(;
     pages=[
         "Home" => "index.md",
         "Zoo" => "zoo.md",
+        "Plotting" => "plotting.md",
         "JuMP" => "jump.md",
         "Gmsh" => "gmsh.md",
+        "PyAMG" => "pyamg.md",
     ],
 )
 

docs/src/index.md

@@ -67,7 +67,8 @@ such as `max_coarse` to the underlying builder):
 - `amg_smoothed_aggregation(; kwargs...)` — smoothed aggregation, via the same package.
 - `amg_pyamg(; solver=:rootnode, kwargs...)` — the Python
   [`pyamg`](https://github.com/pyamg/pyamg) package (`:rootnode` energy-minimization,
-  `:smoothed_aggregation`, or `:ruge_stuben`), imported lazily through `PyCall`.
+  `:smoothed_aggregation`, or `:ruge_stuben`); provided by the
+  [PyAMG extension](pyamg.md) (load PyCall).
 
 ```julia
 mgb_solve(assemble(amg(fem2d_P1(); prolongator = amg_ruge_stuben(max_coarse=4)); p=1.5); verbose=false);
@@ -77,8 +78,7 @@ mgb_solve(assemble(amg(fem2d_P1(); prolongator = amg_ruge_stuben(max_coarse=4));
 
 A 1d p-Laplace problem:
 ```@example 1
-using PyPlot # hide
-using MultiGridBarrier
+using MultiGridBarrier, PyPlot   # PyPlot enables the plotting extension
 nodes = collect(range(-1.0, 1.0, length=33))
 geom = fem1d(; nodes)
 plot(mgb_solve(assemble(amg(geom); p=1.0); verbose=false));
@@ -252,6 +252,12 @@ multigrid hierarchy on the fine mesh. To refine the mesh first, compose with
 The legacy `geometric_mg(geom, L)` builds a geometric-subdivision hierarchy instead of
 AMG; it remains available for callers that specifically want geometric transfers.
 
+Visualization is an opt-in extension: load PyPlot
+(`using MultiGridBarrier, PyPlot`) and `plot` works on every solution and
+geometry — matplotlib for 1D/2D, PyVista for 3D volumes, surfaces, and curves;
+see [Plotting](plotting.md). Without PyPlot the solver core has no Python
+dependency at all.
+
 Prefer stating problems in an algebraic modeling language? The
 [JuMP front end](jump.md) accepts standard `@variable`/`@constraint`/`@objective`
 syntax and lowers it to this same pipeline, building the hierarchy automatically.
@@ -381,7 +387,7 @@ Private = false
 Modules = [MultiGridBarrier]
 Order   = [:type]
 Private = false
-Filter = t -> !(nameof(t) in (:MGBModel, :Coef, :EpiPower, :deriv, :integral, :set_start, :mgb_solution, :solver_log, :On, :Broken, :Continuous, :Uniform, :gmsh_import))
+Filter = t -> !(nameof(t) in (:MGBModel, :Coef, :EpiPower, :deriv, :integral, :set_start, :mgb_solution, :solver_log, :On, :Broken, :Continuous, :Uniform, :gmsh_import, :amg_pyamg, :MGB3DFigure, :HTML5anim))
 ```
 
 # Functions reference
@@ -390,7 +396,7 @@ Filter = t -> !(nameof(t) in (:MGBModel, :Coef, :EpiPower, :deriv, :integral, :s
 Modules = [MultiGridBarrier]
 Order   = [:function]
 Private = false
-Filter = t -> !(nameof(t) in (:MGBModel, :Coef, :EpiPower, :deriv, :integral, :set_start, :mgb_solution, :solver_log, :On, :Broken, :Continuous, :Uniform, :gmsh_import))
+Filter = t -> !(nameof(t) in (:MGBModel, :Coef, :EpiPower, :deriv, :integral, :set_start, :mgb_solution, :solver_log, :On, :Broken, :Continuous, :Uniform, :gmsh_import, :amg_pyamg, :MGB3DFigure, :HTML5anim))
 ```
 
 # Index

docs/src/plotting.md

@@ -0,0 +1,77 @@
+```@meta
+CurrentModule = MultiGridBarrier
+```
+
+# Plotting
+
+Visualization is a package extension (`MultiGridBarrierPyPlotExt`): it loads
+automatically once both MultiGridBarrier and
+[PyPlot](https://github.com/JuliaPy/PyPlot.jl) are imported, and its methods
+extend `PyPlot.plot` and `PyPlot.savefig` — so `plot` is one function whether
+you hand it matplotlib arrays or MultiGridBarrier solutions. Without PyPlot,
+the solver core has no Python dependency at all.
+
+!!! note "Requires PyPlot"
+    Add PyPlot to your environment (`pkg> add PyPlot`) and load both packages:
+    `using MultiGridBarrier, PyPlot`. PyCall arrives with PyPlot; the 3D
+    plotters additionally import the Python `pyvista` package lazily on first
+    use (installing it from conda-forge if necessary).
+
+## What `plot` does
+
+```julia
+plot(sol::MGBSOL, k=1; kwargs...)          # k-th solution component on sol.geometry
+plot(sol::ParabolicSOL, k=1; kwargs...)    # HTML5 animation of a time-dependent solve
+plot(M::Geometry, z::Vector; kwargs...)    # a nodal vector on a geometry
+plot(M::Geometry, ts, U::Matrix; ...)      # animate the columns of U at times ts
+```
+
+The rendering is chosen by the geometry:
+
+- **1D** (`fem1d`, `spectral1d`): line plot. Spectral solutions are sampled
+  through `interpolate`; choose the points with `x = -1:0.01:1`.
+- **2D FEM** (`fem2d_P1`, `fem2d_P2`, `fem2d`): triangulated surface plot
+  (`plot_trisurf`) on the mesh structure.
+- **2D spectral** (`spectral2d`): 3D surface plot on an evaluation grid
+  (`x = -1:0.01:1, y = -1:0.01:1`).
+- **3D FEM** (`fem3d`) and **embedded manifolds** (`fem2d`/`fem1d` with
+  `ambient = Val(3)`): PyVista renders — volume + isosurfaces + slices for
+  volumes, colored quad meshes for surfaces, tubes for curves. These return an
+  [`MGB3DFigure`](@ref) (PNG bytes) that displays inline in Jupyter/Documenter
+  and saves with `savefig(fig, "out.png")`.
+
+Remaining keyword arguments pass through to the underlying PyPlot/PyVista
+calls.
+
+```@example plotting
+using MultiGridBarrier, PyPlot
+sol = mgb_solve(assemble(amg(subdivide(fem2d_P2(), 3)); p = 1.0); verbose = false)
+plot(sol); savefig("plotting_demo.svg"); nothing # hide
+close() # hide
+```
+![](plotting_demo.svg)
+
+## Animations
+
+Time-dependent solutions animate into an [`HTML5anim`](@ref) — an HTML5
+`<video>` that embeds inline in Jupyter, Pluto, and Documenter:
+
+```julia
+sol = parabolic_solve(amg(subdivide(fem2d_P1(), 3)); h = 0.1)
+plot(sol)                       # matplotlib.animation (1D/2D)
+plot(sol.geometry, ts, U)       # or animate any matrix of frames directly
+```
+
+The animation advances at a fixed frame rate (`frame_time` seconds per video
+frame); for irregular `ts`, each video frame shows the latest data frame with
+timestamp ≤ the current video time. In 1D/2D the video is produced by
+`matplotlib.animation` (`embed_limit` caps the embedded size in megabytes, and
+the `printer` callback receives the Matplotlib animation object, e.g. to
+`anim.save("out.mp4")`). In 3D, PyVista frames are encoded by `ffmpeg`.
+
+## Result types
+
+```@docs
+MGB3DFigure
+HTML5anim
+```

docs/src/pyamg.md

@@ -0,0 +1,45 @@
+```@meta
+CurrentModule = MultiGridBarrier
+```
+
+# PyAMG prolongators
+
+The multigrid hierarchy that `amg(geom)` builds is driven by a *prolongator* —
+a callable mapping a stiffness matrix to its level prolongations. Two
+pure-Julia factories ship in the core (`amg_ruge_stuben`, the default, and
+`amg_smoothed_aggregation`, both via
+[AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl)).
+A third is backed by the Python [pyamg](https://github.com/pyamg/pyamg)
+package and lives in the `MultiGridBarrierPyAMGExt` extension:
+
+!!! note "Requires PyCall"
+    Add PyCall to your environment (`pkg> add PyCall`) and load both packages:
+    `using MultiGridBarrier, PyCall`. The Python `pyamg` and `scipy` packages
+    are imported lazily on the first call, installing from conda-forge if
+    necessary. (Loading PyPlot for [plotting](plotting.md) also loads PyCall,
+    so the plotting setup enables this extension too.)
+
+```julia
+using MultiGridBarrier, PyCall
+geom = subdivide(fem2d_P1(), 4)
+mg   = amg(geom; prolongator = amg_pyamg(solver = :rootnode))
+sol  = mgb_solve(assemble(mg; p = 1.0))
+```
+
+`solver` selects pyamg's `:rootnode` (energy-minimization, the default),
+`:smoothed_aggregation`, or `:ruge_stuben`; remaining keyword arguments are
+forwarded to the pyamg solver constructor.
+
+When is this worth reaching for? The pure-Julia Ruge–Stüben default is robust
+across the package's benchmarks. Aggregation-based coarsening —
+`amg_pyamg(solver = :rootnode)` in particular — is an escape hatch for highly
+anisotropic problems near `p = 1`, especially combined with the
+`auxiliary_postprocess` keyword of `amg` (see the [`amg`](@ref) docstring),
+where classical coarsening can blow up Newton iteration counts on the central
+path.
+
+## API reference
+
+```@docs
+amg_pyamg
+```

docs/src/zoo.md

@@ -20,8 +20,7 @@ dimension, on the coarsest mesh that still draws an interesting picture.
 Tolerances are loose (`tol = 1e-3`) and progress bars are suppressed.
 
 ```@example zoo
-using PyPlot  # hide
-using MultiGridBarrier
+using MultiGridBarrier, PyPlot   # PyPlot enables the plotting extension
 nodes1 = collect(range(-1.0, 1.0, length=17))     # 1D
 mg1 = amg(fem1d(; nodes=nodes1))
 mg2 = amg(subdivide(fem2d_P1(), 3))               # 2D

ext/MultiGridBarrierPyAMGExt.jl

@@ -0,0 +1,51 @@
+# MultiGridBarrierPyAMGExt — the pyamg-backed AMG prolongator, as a package
+# extension. Loads automatically when both MultiGridBarrier and PyCall are
+# imported (`using MultiGridBarrier, PyCall`); supplies the method of the
+# parent stub `amg_pyamg` (src/amg_prolongators.jl). The Python pyamg and scipy
+# packages are imported lazily on the first call, installing from conda-forge
+# if necessary.
+module MultiGridBarrierPyAMGExt
+
+using PyCall
+using SparseArrays
+import MultiGridBarrier
+import MultiGridBarrier: amg_pyamg
+
+# Convert a scipy sparse matrix (PyObject) to a Julia SparseMatrixCSC{Float64,Int}.
+function _scipy_to_julia(Pyobj)
+    csc     = Pyobj.tocsc()
+    indptr  = Vector{Int}(csc.indptr)      # 0-based CSC column pointers
+    indices = Vector{Int}(csc.indices)     # 0-based row indices
+    data    = Vector{Float64}(csc.data)
+    shape   = csc.shape
+    m, n    = Int(shape[1]), Int(shape[2])
+    colptr  = indptr .+ 1
+    rowval  = indices .+ 1
+    return SparseMatrixCSC{Float64,Int}(m, n, colptr, rowval, data)
+end
+
+function MultiGridBarrier.amg_pyamg(; solver::Symbol=:rootnode, kwargs...)
+    return function (K64::SparseMatrixCSC{Float64,Int})
+        pyamg = pyimport_conda("pyamg", "pyamg", "conda-forge")
+        scipy_sparse = pyimport_conda("scipy.sparse", "scipy", "conda-forge")
+        # Build a scipy CSC matrix from the Julia CSC arrays (0-based indices).
+        A = scipy_sparse.csc_matrix(
+            (K64.nzval, K64.rowval .- 1, K64.colptr .- 1),
+            shape = (size(K64, 1), size(K64, 2)))
+        ml = if solver === :rootnode
+            pyamg.rootnode_solver(A; kwargs...)
+        elseif solver === :smoothed_aggregation
+            pyamg.smoothed_aggregation_solver(A; kwargs...)
+        elseif solver === :ruge_stuben
+            pyamg.ruge_stuben_solver(A; kwargs...)
+        else
+            throw(ArgumentError("amg_pyamg: unknown solver $(solver); " *
+                "use :rootnode, :smoothed_aggregation, or :ruge_stuben"))
+        end
+        # The coarsest level carries no P; take all levels but the last.
+        levels = collect(ml.levels)
+        return [_scipy_to_julia(lvl.P) for lvl in levels[1:end-1]]
+    end
+end
+
+end # module