BpjSrSSKrSSKJr SSKJrJrJrJ r SSK J r SSK J r SSKJr S!S jrS r"S S 5rS rS"Sjr"SS5rSSSS.SjrSrSrS#SjrSrSrSSSSSSSS.SjrSSSSSS.S jrg)$u Regrid (2-D smoothing B-splines via separable 1-D FITPACK kernels) ================================================================== 1) Overview ----------- This module fits a bivariate tensor-product B-spline surface to gridded data `Z[i, j]` sampled on strictly increasing coordinates `x[i]`, `y[j]`. It mirrors the adaptive spirit of FITPACK's REGRID/fpgrre (knot growth + smoothing parameter search to meet a target residual `s`) while keeping the control flow explicit in Python and returning a 2-D `NdBSpline`. **Key conventions** - **Penalty scaling is 1/p** - *same as FITPACK*. The augmented system stacks `(D / p)` under the data matrix. Smaller `p` -> stronger smoothing; larger `p` -> weaker smoothing (approaching interpolation). - `p == -1` is used **as a sentinel for `p = ∞`** (interpolatory limit): when seen, the solver omits penalty rows entirely. 2) Mathematics (and how this differs from FITPACK's REGRID) ----------------------------------------------------------- Let `A_x`, `A_y` be banded 1-D design matrices and `D_x`, `D_y` the banded 1-D roughness (difference) matrices returned by FITPACK/Dierckx's 1-D APIs (`data_matrix`, `disc`). The 2-D smoothing objective is minimize ||A c - z||^2 + (1/p) ||D c||^2, (1) implemented by the **augmented** least squares system [ A ] c ~ [ z ] [D/p] [ 0 ]. **Same as FITPACK:** Equation (1) uses the *1/p* convention. **Different in this module:** We realize the 2-D problem by composing **1-D** banded operators in two separable passes (x then y), instead of calling a monolithic 2-D routine. This makes the mathematics transparent while producing the same normal equations structure that REGRID targets internally. 3) Residual energy: definition and use -------------------------------------- After solving on the current knots (initially at the interpolatory limit, `p = ∞`), we evaluate the surface and form R = Z - Zhat, fp = sum(R[i, j]^2). We then **project** residual energy to knot spans along each axis: - Row energy `row_energy[i] = sum(R[i, j]^2, j)` -> accumulated into `fpintx`. - Column energy `col_energy[j] = sum(R[i, j]^2, i)` -> accumulated into `fpinty`. These per-span energies guide **adaptive knot insertion**: the algorithm picks high-energy spans (skipping zero-length ones) and inserts data-aligned knots (e.g., median sample within the span), with simple batch sizing and headroom limits. This is conceptually the same idea as REGRID's `fpint` arrays; here it is implemented explicitly and vectorized in Python for clarity. 4) Solver used here vs. FITPACK's REGRID ---------------------------------------- **This module (separable 1-D composition):** - Uses **1-D FITPACK/Dierckx kernels** (`data_matrix`, `disc`, `qr_reduce`, `fpback`) to solve 2-D via two passes: 1) augment/QR/backsolve along **x**, producing an intermediate; 2) augment/QR/backsolve along **y**, producing the coefficient grid `C`. - The augmented rows are stacked as `[A; D/p]` (1/p scaling, **same as FITPACK**). - If `fp > s` after knot growth, performs a **scalar search in `p`** using a ratio-of-roots routine; `p == -1` is treated as `p = ∞` for the interpolatory reference without penalty rows. **FITPACK REGRID (monolithic 2-D routine):** - Implements the same mathematical objective and **1/p** penalty scaling, but inside a specialized, fully 2-D Fortran routine with in-place Givens/QR and a rational update for `p`. - Handles residual partitioning, knot insertion, and `p` updates internally. **Practical difference:** We build the 2-D solve from **1-D building blocks**, which keeps each stage observable/testable and uses the same low-level kernels as FITPACK, but without relying on the single monolithic REGRID entry point. 5) Execution flow (who calls whom) ---------------------------------- **Top-level** 1. `_regrid` - Validates shapes/monotonicity and normalizes `bbox`. - Dispatches to `_regrid_fitpack(...)`. **Core driver** 2. `_regrid_fitpack` - Clips inputs to `bbox` via `_apply_bbox_grid` -> `(x_fit, y_fit, Z_fit)`. - If `s == 0` (interpolation path): - Build initial not-a-knot vectors: `tx = _not_a_knot(x_fit, kx)`, `ty = _not_a_knot(y_fit, ky)`. - Build design matrices: `(Ax, Ay, Q) = build_design_matrices(...)`. - Solve once at `p = -1` (inf): `C0, fp = _solve_2d_fitpack(...)`. - Return `return_NdBSpline(...)`. - Else (`s > 0`, smoothing with adaptive knots): - Initialize clamped no-interior-knot vectors with `_initialise_knots` (for x and y). - **Knot-growth loop** (up to `len(x_fit)+len(y_fit)` iterations): 1) Build design matrices: `(Ax, Ay, Q) = build_design_matrices(...)`. 2) Interpolatory reference solve (`p = -1`): `C0, fp = _solve_2d_fitpack(...)`. 3) If both `tx`, `ty` are at minimal size, record `fp0 = fp` (used to bracket the p-search). 4) If `fp < s`: **stop growth** and proceed to p-search / finalize. 5) Compute residuals for knot insertion: - Convert packed to CSR once per axis for evaluation: `_Ax = Ax.tocsr(...)`, `_Ay = Ay.tocsr(...)` - Evaluate `Z0 = _Ax @ C0 @ _Ay.T`, residual `R = Z_fit - Z0`. 6) Insert knots on alternating axes using energy heuristics: - If last axis was `y`, grow `tx` with `_add_knots(..., residuals=np.sum(R**2, axis=1), ...)`. - Else, grow `ty` with `_add_knots(..., residuals=np.sum(R**2, axis=0), ...)`. - Update `fpold`, `nplus{ x|y }`, `last_axis`. - If growth ended with both axes still minimal: return `return_NdBSpline(...)`. - **Finite-p smoothing (if needed):** - Build 1-D penalty operators: `(Drx, _, ncx) = disc(tx, kx)`, `(Dry, _, ncy) = disc(ty, ky)`, wrap as `PackedMatrix`. - Rebuild design matrices on the final `(tx, ty)`. - Run `_p_search_hit_s(...)` to find `p*` such that `fp(p*) ~ s`: - Internally constructs `F(...)` which evaluates `fp(p)` by calling `_solve_2d_fitpack(...)`. - Uses `root_rati` on `g(p) = fp(p) - s` with the interpolatory reference at `p = inf` and `fp0` bracket. - Return `return_NdBSpline(...)`. **Separable solver (used by both interpolation and p-search)** 3. `_solve_2d_fitpack` - Forms augmented banded systems via `_stack_augmented_fitpack`: - X-side: `Ax_aug = [Ax; Dx/p]` if `p != -1`, else just `Ax`. - Y-side: `Ay_aug = [Ay; Dy/p]` if `p != -1`, else just `Ay`. - Pads RHS `Q` with zeros when penalties are stacked. - **Stage X**: `_dierckx.qr_reduce(Ax_aug, ...)` then `_dierckx.fpback(...)` -> intermediate `T`. - Transpose `T` to feed Y solves column-wise; pad if needed. - **Stage Y**: `_dierckx.qr_reduce(Ay_aug, ...)` then `_dierckx.fpback(...)` -> coefficients `C`. - Build CSR design matrices once: `_Ax = Ax.tocsr(...)`, `_Ay = Ay.tocsr(...)`. - Evaluate `zhat = _Ax @ C.T @ _Ay.T`, compute `fp = ||Z_fit - zhat||^2`. - Return `C.T` (in `(nx_coef, ny_coef)` layout) and `fp`. **Helpers** - `_apply_bbox_grid(...)` - slices to `(x_fit, y_fit, Z_fit)`. - `build_design_matrices(...)` - wraps `_dierckx.data_matrix` and returns `PackedMatrix` wrappers + `Q`. - `_initialise_knots(...)`, `_add_knots(...)` - FITPACK-style knot bookkeeping/growth. - `disc(...)` - 1-D roughness operators (packed band form). - `F` + `root_rati` - maps `p -> fp(p)` and finds `p*` with `fp(p*) ~ s`. - `return_NdBSpline(...)` - packs `(tx, ty, C)` into an `NdBSpline`. N) NdBSpline) root_ratidiscadd_knot _not_a_knot)_dierckx) csr_array) _validate_intcz[UR5S:wa [S5e[US5n[US5nUS:dUS:a [S5eURR SSn[ R"U5n[ R"U5nU(Ga@URS:XdURS:XaD[ R"URUR4U-URRS9nU$URS:a=[ R"[ R"U5S :5(d [S 5eURS:a=[ R"[ R"U5S :5(d [S 5e[ R"XS S 9up[ R"X4SS9n U"XU4URS9nU$UR UR :wa[ R "X5upURS:Xa6[ R"UR U-URRS9$[ R"UR#5UR#54SS9n U"XU4URS9nUR%UR U-5$)a Evaluate a 2D `NdBSpline` like a classical bivariate API. Parameters ---------- ndbs : NdBSpline A 2D spline object (``len(ndbs.t) == 2``). x, y : array_like Sample locations. If ``grid=True``, these must be 1-D strictly increasing vectors. If ``grid=False``, they can be broadcastable arrays of the same shape. dx, dy : int, optional Derivative orders along `x` and `y` respectively, by default 0. grid : bool, optional If True, evaluate on the cartesian product of `x` and `y`; otherwise treat `(x, y)` as paired coordinates, by default True. Returns ------- ndarray or (ndarray, dict) Evaluated values with shape: - ``(len(x), len(y), ...)`` if ``grid=True``. - ``x.shape + ...`` if ``grid=False``. Raises ------ ValueError If `ndbs` is not 2D, derivatives are negative, or monotonicity checks fail. Notes ----- This is a thin convenience wrapper around ``NdBSpline.__call__`` with input validation and optional profiling. z*ndbs must be a 2D NdBSpline (len(t) == 2).dxdyrz,order of derivative must be positive or zeroNdtypez1x must be strictly increasing when `grid` is Truez1y must be strictly increasing when `grid` is Trueij)indexingaxis)nu extrapolate)lent ValueErrorr cshapenpasarraysizezerosralldiffmeshgridstackrbroadcast_arraysravelreshape) ndbsxyrrgridtrailingvalsXYxis U/var/www/html/pdf-tiff/venv/lib/python3.13/site-packages/scipy/interpolate/_regrid.py_ndbspline_call_like_bivariater4sF 466{aEFF r4 B r4 B AvaGHHvv||ABH 1 A 1 A 66Q;!&&A+88QVVQVV,x7tvv||LDK FFaK"&&s):";";PQ Q FFaK"&&s):";";PQ Q{{1$/ XXqf2 &B81A1AB 77agg &&q,DA 66Q;88AGGh.dffllC C XXqwwy!''),2 6B81A1AB||AGGh.//c[US5[US5pCUupVUSRX5- S- XF- S- 5n[USUS4Xr5$)aO Build a 2D ``NdBSpline`` from knot vectors and a coefficient grid. Parameters ---------- fp : float Residual sum of squares of the produced fit (kept for upstream use). tck : tuple Tuple ``(tx, ty, C)`` where ``tx``, ``ty`` are knot vectors and ``C`` is a coefficient array with shape ``(nx - kx - 1, ny - ky - 1)`` or a compatible shape that can be reshaped to that. degrees : tuple of int Degrees ``(kx, ky)`` along x and y. Returns ------- NdBSpline The constructed 2D spline. Notes ----- Only repacks the coefficient grid; ``fp`` is not used internally here. rrr )rr)r)fptckdegreesnxnykxkyrs r3return_NdBSpliner>s\0Q[#c!f+ FB Arw{BGaK0A c!fc!f%q 22r5c:\rSrSrSrSr\S5rSrSr Sr g) PackedMatrixia_A simplified CSR format for when non-zeros in each row are consecutive. Assuming that each row of an `(m, nc)` matrix 1) only has `nz` non-zeros, and 2) these non-zeros are consecutive, we only store an `(m, nz)` matrix of non-zeros and a 1D array of row offsets. This way, a row `i` of the original matrix A is ``A[i, offset[i]: offset[i] + nz]``. c(XlX lX0lgN)aoffsetnc)selfrCrDrEs r3__init__PackedMatrix.__init__s r5cLURRSUR4$)Nr)rCrrE)rFs r3rPackedMatrix.shapesvv||A''r5cn[R"UR5nURRSn[ URS5H_n[ UR URU- U5nURUSU24XURUURUU-24'Ma U$)Nrr)rr"rrCrangeminrErD)rFoutneleminels r3todensePackedMatrix.todenseshhtzz" Qsyy|$Adgg A.6C:>&&DSD/C4;;q>$++a.3"666 7% r5cURR5n[R"URUS-5R SUS-5nU[R "US-URRS9-nUR5n[R "SUS-US--US-URRS9n[XEU4X#U- S- 4S9$)Nrrrr)r) rCr(rrepeatrDr)arangerr )rFkmlen_tdataindicesindptrs r3tocsrPackedMatrix.tocsr svv||~))DKK1-55b!A#>BIIac1B1BCC--/1q1uQ/Q!%!2!24 F #ai!m$ r5)rCrErDN) __name__ __module__ __qualname____firstlineno____doc__rGpropertyrrRr]__static_attributes__r5r3r@r@s* ((r5r@cUS:Xa5URR5URR5U4$US-n[R"X!R S-US-4[ S9nURSU2SS24USU2SU24'URU- XbS2SS24'[R"URUR45nXgU4$)a Builds augmented banded matrix. Parameters ---------- A : PackedMatrix Banded data/design matrix for one axis (from `_dierckx.data_matrix`). D : PackedMatrix Banded roughness (difference) penalty matrix for the same axis (from `disc`). nc : int Number of top (data) rows from `A` to include. k : int Spline degree (used only for sizing in the current implementation). p : float Smoothing parameter. The effective penalization term is scaled as **1/p**: larger `p` means *less* smoothing (approaching interpolation). If `p == -1`, it signals *p -> inf*, i.e. a pure interpolatory system with **no** penalty rows appended. Returns ------- AA : ndarray Augmented banded matrix with `A` stacked over `(D / p)` when `p != -1`. offset : ndarray Concatenated band offsets for the augmented matrix. nc : int Returned unchanged for downstream convenience. rrrr rN)rCcopyrDrr"rfloat concatenate)ADrErWpnzAArDs r3_stack_augmented_fitpackrp2s< Bwssxxz188==?B.. QB 2 ?AE*% 8B33ssAv;BssCRCxLqBsAvJ ^^QXXqxx0 1F r>r5c [R"U5n [R"U 5n[X URSXC5unnnURn[XURSXs5unnnURnUS:waK[R "U[R "U RSURS4[S9/5n[R"UUUU5 [R"UUUX%XMUS5 un n[R"UR5nUS:waK[R "U[R "U RSURS4[S9/5n[R"UUUU5 [R"UUXXUUS5 unnnURXFRS[U55nURXyRS[U55nUUR-UR-nU U- n[R"[R "U55nURUU4$)a Solve the 2-D tensor-product spline system using separable banded QR. ================================================================ Mathematical model (step by step, plain text) ================================================================ Shapes: Z : (mx, my) -> original data Ax, Ay : design matrices for x and y Dx, Dy : roughness penalty matrices for x and y C : (nx, ny) -> spline coefficients to solve for Surface approximation: Zhat = Ax * C * Ay^T Objective (smoothing formulation): minimize ||Ax*C*Ay^T - Z||^2 + (1/p)*(||Dx*C||^2 + ||C*Dy^T||^2) In practice (FITPACK-style separable approach), we solve this in two stages: -------------------------------------------------------- Stage 1 (x-direction solve for all y-columns together): -------------------------------------------------------- For each column of Z: minimize ||Ax*T - Z||^2 + (1/p)*||Dx*T||^2 This is equivalent to the augmented least-squares system: [Ax] [Z] [Dx/p] * T = [0] i.e. minimize || [Ax; Dx/p]*T - [Z; 0] ||^2 The solution T is obtained by QR reduction and back-substitution. -------------------------------------------------------- Stage 2 (y-direction solve using transposed result): -------------------------------------------------------- Now treat T^T as the new RHS for the y-direction: minimize ||Ay*C^T - T^T||^2 + (1/p)*||Dy*C^T||^2 Equivalent to augmented system: [Ay] [T^T] [Dy/p] * C^T = [0] i.e. minimize || [Ay; Dy/p]*C^T - [T^T; 0] ||^2 Solving this gives C^T (then transposed back to C). -------------------------------------------------------- Interpolation limit: -------------------------------------------------------- If p == -1, penalties are omitted (Dx, Dy are not stacked). The solver behaves as a near-interpolating system. -------------------------------------------------------- Residual computation: -------------------------------------------------------- Zhat = Ax * C * Ay^T R = Z - Zhat fp = sum(R^2) Parameters ---------- Ax, Ay : PackedMatrix Banded data matrices for the x and y axes. Q : ndarray, shape (mx, my) RHS data grid (copied from `Z`). p : float Smoothing parameter. The penalty term is scaled as **1/p**. Setting `p == -1` signals *p -> inf* (interpolation, omit penalty). kx, ky : int Spline degrees along x and y. tx, ty : ndarray Knot vectors along x and y. x_x, x_y : ndarray Sample coordinates. z : ndarray Original data grid for residual evaluation. Dx, Dy : ndarray Banded roughness penalty matrices for x and y. Optional, Only needed when ``p != -1``. Returns ------- C : ndarray 2-D B-spline coefficient grid. fp : float Residual sum of squares between fitted surface and `z`. R : ndarray, shape (mx, my) Residual matrix ``z - zhat``, where ``zhat = Ax @ C @ Ay.T``. Notes ----- This performs two separable QR solves (x then y), each augmented by `(D / p)` when `p != -1`. Setting `p = -1` skips all penalty rows, yielding an interpolatory surface. The resulting coefficients and residual follow the same conventions as FITPACK's `fpgrre`. rrrrF)r ones_likerprrEvstackr"rir qr_reducefpbackascontiguousarrayTr]rsumsquare)AxAyQrmr<txx_xr=tyx_yzDxDyw_xw_yAx_aug offset_aug_xnc_augxnc_xAy_aug offset_aug_ync_augync_yrw_Cr7_Ax_AyzhatRs r3_solve_2d_fitpackrZsR ,,s C ,,s C%=  R%$!FL' 55D%=  R%$!FL' 55D Bw IIq"((BHHQK#T"U5T- $rBrf)rmfp_atss r3g_p_search_hit_s..gsQx!|r5rrg-q=)r)rrinfmaxrrootr)rzrr{rr|r<r}r~r=rrrrfp0rrrrfpmsbracketftolrp_starfp_starC_starrs ` @r3_p_search_hit_srhsl bbars 'E R5DS1W~~.G q7{E "D!VWd%8A VVFFmG WWF 67 ""r5c[UVs/sHoDSLPM sn5(aXU[S5[S54$UupVpxXV:aXx:d [S5e[R"X:X:*-5Sn [R"X:X:*-5Sn U R S:XdU R S:Xa [S5eX XU[R "X5[RU [RU 4$s snf)a Restrict (x, y, Z) to a rectangular bounding box. Parameters ---------- x, y : ndarray Monotonic sample coordinates. Z : ndarray, shape (len(x), len(y)) Data grid. bbox : sequence of 4 scalars or None ``(xb, xe, yb, ye)``; any element may be None to skip clipping. Returns ------- x_fit, y_fit, Z_fit : ndarray Sliced arrays restricted to bbox. ix, iy : slice or ndarray Indexers mapping from full arrays to the restricted ones. Raises ------ ValueError If bbox is invalid or excludes all samples along an axis. Nz%bbox must satisfy xb < xe and yb < yerz$bbox excludes all samples in x or y.)r#slicerrwherer!ix_s_) r+r,Zbboxbboxixbxeybyeixiys r3_apply_bbox_gridrs2 t ,teTMt ,--Qd U4[00NBB G@AA 17qw' ( +B 17qw' ( +B ww!|rww!|?@@ 5!%266">*BEE"IruuRy @@ -s Dc[R"U5n[R"U5n[R"XXW5upn [R"XXh5upnUR 5n[ XU 5[ XU5U4$rB)rrrr data_matrixrhr@)r+r,rr}rr<r=rrrzoffset_xrr{offset_yrr|s r3_build_design_matricesrsv ,,q/C ,,q/C!--aR=B$!--aR=B$ A t , t ,  r5cUc[X-S-SU-S-5nO%USUS--:a[SU<SSUS--S35eSUS--nX-S-n[R"U/US--U/US---5nXtXV4$)a Initialize a non-periodic knot vector. Parameters ---------- m : int Number of data points (equivalent to len(x) if x were provided). xb, xe : float Domain endpoints used to seed the initial knot vector with no internal knots. k : int Spline degree. nest : int, optional Storage cap for knots. If None, defaults to max(m + k + 1, 2*k + 3). Must satisfy nest >= 2*(k + 1); otherwise a ValueError is raised. Returns ------- t : 1-D ndarray Initial knot vector with no internal knots: [xb]*(k+1) + [xe]*(k+1). nest : int The finalized storage cap for knots. nmin : int Lower bound on knot count. nmax : int Upper bound on knot count. What this does -------------- - Computes defaults and bounds used by FITPACK-style knot growth: * nest: storage cap for knots (defaults to max(m + k + 1, 2*k + 3)) * nmin: minimal knot count (2*(k+1)) * nmax: maximal knot count (m + k + 1) - Returns an initial knot vector with no internal knots: t = [xb]*(k+1) + [xe]*(k+1) rr z`nest` too small: nest = z < 2*(k+1) = .)rrrr )rXrrrWnestnminnmaxrs r3_initialise_knotsrsH |1519acAg& !QU) 9$-1Q3yPQRS S a!e9D 519D B41:ac *+A D r5c TU[-n URn Xr- n X:XaSn O:X- nX:a[X-U- 5OU S-n[U S-[ XS-S55n [ U 5H<n[ XX5nURSn X:a[X5U 4s $X:dM9X:4s $ X:4$)aE Knot-growth helper for knot-finding loop (non-periodic). Parameters ---------- x : 1-D ndarray Strictly increasing sample coordinates. k : int Spline degree. s : float Target smoothing. t : 1-D ndarray Current knot vector to be grown. nmin, nmax : int Lower/upper bounds on knot count (from initialisation). nest : int Storage cap for total knots. fp, fpold : float Current and previous residual sums of squares. Used to update nplus. residuals : 1-D ndarray Most recent residual signal used by `add_knot` to decide placement. nplus : int Previous iteration's proposed number of knots; used to update the next nplus. Returns ------- t_new, nplus : tuple Updated knot vector and the nplus chosen for this step. If n >= nmax, t_new is a not-a-knot layout. If n >= nest, t_new is the current vector respecting the storage cap. What this function does ----------------------- - Assumes the caller has already decided to GROW (i.e., checks like |fp - s| < acc or fp < s has FAILED). - Updates nplus (how many knots to add next) using the FITPACK heuristic based on the previous improvement (delta = fpold - fp). - Inserts up to nplus new internal knots using `add_knot(x, t, k, residuals)`. - Stops early if storage or interpolation caps are reached: * if n >= nmax: switch to interpolation layout (not-a-knot) and return * if n >= nest: return current t respecting storage cap How it compares with _fitpack_repro.py::_generate_knots_impl ------------------------------------------------------------- Similarities: 1) Same growth logic for nplus: - Use delta = fpold - fp with ratio fpms/delta - Apply min/max caps (doubling and halving behavior) 3) Same storage guard: - If n reaches nest, stop and return current t 4) Same end behavior at the "interpolating" cap: - When n >= nmax, switch to not-a-knot layout and return Differences: 1) API style: - _generate_knots_impl is a generator that yields trial knot vectors and recomputes residuals/fp internally on each iteration. - `_add_knots` is a stateful helper that only grows knots; it expects the caller to handle residual computation and fp/fpold updates between calls. 2) Periodicity: - _generate_knots_impl supports periodic=True. - `_add_knots` is non-periodic only; it uses not-a-knot when n >= nmax. 3) Residual computation: - _generate_knots_impl calls an internal residual routine each iteration. - `_add_knots` does not compute residuals; the caller must supply: residuals (used by add_knot), fp, fpold. 4) Return values: - _generate_knots_impl yields multiple t's and eventually returns None. - `_add_knots` returns: * (t_new, nplus) after inserting knots, * (not_a_knot_t, nplus) if n >= nmax, * (t, nplus) if n >= nest (storage cap). rr r) TOLr!intrMrrLrrr)r+rWrrrrrr7fpold residualsnplusaccnrdeltanpl1js r3 _add_knotsrsZ c'C A 6D  y ,1Ks5<%'(U1WE!GSax34 5\ Q1 ( GGAJ 9q$e+ + 98O'* 8Or5rr2r<r=rrnestxnestyrcL[XX)5upn n U RUS-:dU RUS-:a8[SUSUSUS-SUS-SU RSU RS3 5e[U ScU SOU S5n[U ScU S OU S5n[U S cU SOU S 5n[U S cU S OU S 5nS nUS :XabUcUb [S5e[ X5n[ X5n[ XUUUX45unnn[ UUUUUUXUX5 unnn [UUUU4X445$[U RXX7S9unnnn[U RUUXHS9unnnnS nSn[U5[U5-n S n!S n"S n#[U 5Hn [ XUUUX45unnn[ UUUUUUU UUU U 5 unnn$[U5U:Xa[U5U:XaUn!UU:a OUS:Xa.[XUUUUUUU[R"U$S -SS9U"S9 unn"SnO-[XUUUUUUU[R"U$S -SS9U#S9 unn#Sn[U5U:a[U5U:a OUnM [U5U:Xa [U5U:Xa[WUUW4X445$Sn[UU5un%n&n'[UU5un(n)n*[U%U&U'5n%[U(U)U*5n([ XUUUX45unnn[!UU%UU(UUUXUXUU!UUS9un n+n,[U,UUU+4X445$)a Core adaptive bivariate spline fitter using the 1/p-penalty convention. Parameters ---------- x, y : array_like Strictly increasing coordinate vectors. Z : array_like, shape (len(x), len(y)) Data grid. kx, ky : int, optional Spline degrees along x and y, default 3 (cubic). s : float, optional Target residual (`fp` target). `s = 0` requests an interpolatory surface; `s > 0` triggers smoothing with penalty weight **1/p**. maxit : int, optional Maximum iterations for the `p`-search when smoothing, default 50. nestx, nesty : int or None Max coefficient counts per axis (nesting limits). bbox : sequence of 4 scalars Optional domain limits `(xb, xe, yb, ye)`. Use `None` entries to skip. Returns ------- NdBSpline Fitted 2-D spline surface. Notes ----- The internal smoothing parameter `p` follows the **inverse**-penalty rule: penalty term is 1/p. Hence, larger `p` -> weaker smoothing (approaching interpolation), while smaller `p` -> stronger smoothing. A sentinel value `p == -1` is interpreted as *p = inf*, corresponding to an exact (interpolatory) fit. The iterative process adaptively grows knot vectors based on residual energy and optionally performs a 1-D search over `p` to satisfy `fp ~ s`. rz/Not enough samples inside bbox for degrees (kx=z, ky=z ). Need at least k+1 per axis: (z, z). Got (z).rNrr rrzs == 0 is interpolation only)rr,r)rrrr7rrrr+)rr)rr!rrirrrr>rrrLrrrxrr@r)-r+r,rr<r=rrrrrx_fity_fitZ_fitrrrrrrmr}rrzr{r|C0r7nminxnmaxxnminynmaxyr last_axismpmrnplusxnplusyrDrx offset_dxnc_dxDry offset_dync_dyC_smfp_sms- r3_regrid_fitpackrsR!1q ?E%A zzR!V b1f 5=bTrdK,,.qD6BqD6:JJB8 <= 0``. The underlying formulation uses a **1/p** penalty, meaning: - small `p` -> heavy smoothing, - large `p` -> light smoothing (approaching interpolation). Setting `p == -1` internally denotes *p = inf*, i.e. a pure interpolant. maxit : int, optional Maximum iterations for `p`-search if invoked. Returns ------- NdBSpline Fitted bivariate spline surface. Nrrzx must be strictly increasingzy must be strictly increasingrz7x dimension of z must have same number of elements as xrz7y dimension of z must have same number of elements as yzs should be s >= 0.0)rz"bbox shape should be (4,), found: r) rr rir(r#r$rr!rr)r+r,rrr<r=rrs r3_regridr!s=8 |vax 1E"A 1E"A 1E"A 88D>D aA 66"''!*s" # #899 66"''!*s" # #899vvRSSvvRSS C/00 :: =djj\JKK  a2%$T ++r5)rrT)NNrB)rcnumpyrscipy.interpolate._ndbspliner scipy.interpolate._fitpack_reprorrrrr scipy.sparser scipy._lib._utilr r4r>r@rprrrrrrrrrrrfr5r3rsTj2,,"*G0T3<))X&V#'JX??J BF#R%AP 0f tpc D K=\"aAB4+r5