+ nDjRt^RIt^RIHt^RIHtHtHtH t ^RI H t ^RI H t ^RIHtR"RltRt!R R 4tR tR#R lt!R R4tRRRRR^(/RltRtRtR$RltRtRtR^R^RRR^2RRRRRR/R ltRRR^R^RRR^2/R!ltR#)%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_intc\VP4^8wd \R4h\VR4p\VR4pV^8gV^8d \R4hVPP R,p\ P!V4p\ P!V4pV'EdWVP^8XgVP^8XdL\ P!VPVP3V,VPPR7pV#VP^8d?\ P!\ P!V4R84'g \R4hVP^8d?\ P!\ P!V4R84'g \R 4h\ P!WR R 7wr\ P!W3RR 7p V!WV3VPR 7pV#VP VP 8wd\ P !W4wrVP^8Xd>\ P!VP V,VPPR7#\ P!VP#4VP#43RR 7p V!WV3VPR 7pVP%VP V,4#)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).dxdyz,order of derivative must be positive or zero:NNdtypez1x must be strictly increasing when `grid` is Truez1y must be strictly increasing when `grid` is Trueij)indexingaxis)nu extrapolate)lent ValueErrorr cshapenpasarraysizezerosralldiffmeshgridstackrbroadcast_arraysravelreshape) ndbsxyr r gridtrailingvalsXYxis &&&&&& S/data/cameron/venvs/s3viz/lib/python3.14/site-packages/scipy/interpolate/_regrid.py_ndbspline_call_like_bivariater2sF 466{aEFF r4 B r4 B AvaGHHvv||BH 1 A 1 A t 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\V^,4\V^,4rCVwrVV^,PW5, ^, WF, ^, 4p\V^,V^,3Wr4#)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. )rr'r)fptckdegreesnxnykxkyrs&&& r1return_NdBSpliner<s\0Q[#c!f+ FB Arw{BGaK0A c!fc!f%q 22r3cFa]tRtRtoRtRt]R4tRtRt Rt Vt R#) 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*WnW nW0nR#N)aoffsetnc)selfrArBrCs&&&&r1__init__PackedMatrix.__init__s r3cTVPP^,VP3#))rArrC)rDs&r1rPackedMatrix.shapesvv||A''r3c\P!VP4pVPP^,p\ VP^,4Fwp\ VP VPV,, V4pVPVRV13,WVPV,VPV,V,13&Ky V#)N)rr rrArangeminrCrB)rDoutneleminels& r1todensePackedMatrix.todenseshhtzz" Qsyy|$Adgg A.6C:>&&DSD/C4;;q>$++a.3"666 7% r3cVPP4p\P!VPV^,4P RV^,4pV\P !V^,VPPR7,pVP4p\P !^V^,V^,,V^,VPPR7p\WEV3W#V, ^, 3R7#)rKr)rr) rAr&rrepeatrBr'arangerr)rDkmlen_tdataindicesindptrs&&&& r1tocsrPackedMatrix.tocsr svv||~))DKK1-55b!A#>BIIac1B1BCC--/1q1uQ/Q!%!2!24 F #ai!m$ r3)rArCrBN) __name__ __module__ __qualname____firstlineno____doc__rEpropertyrrRr]__static_attributes____classdictcell__ __classdict__s@r1r>r>s2 ((r3r>cVR8Xd6VPP4VPP4V3#V^,p\P!W!P ^,,V^,3\ R7pVPRV1R3,VRV1RV13&VPV, WbR1R3&\P!VPVP34pWgV3#)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. rN:NNNr)rAcopyrBrr rfloat concatenate)ADrCrWpnzAArBs&&&&& r1_stack_augmented_fitpackrr2s< Bwssxxz188==?B.. QB 2 ?AE*% 8B33ssAv;BssCRCxLqBsAvJ ^^QXXqxx0 1F r>r3c &\P!V4p \P!V 4p\W VP^,WC4wpppVPp\WVP^,Ws4wpppVPpVR8wdW\P !V\P !V P^,VP^,3\R7.4p\P!VVVV4\P!VVVW%WMVR4 wpp\P!VP4pVR8wdW\P !V\P !V P^,VP^,3\R7.4p\P!VVVV4\P!VVWWVVR4 wpppVPWFP^,\V44pVPWyP^,\V44pVVP,VP,pV V, p\P!\P !V44pVPVV3#)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`. rFr)r ones_likerrrrCvstackr rkr qr_reducefpbackascontiguousarrayTr]rsumsquare)AxAyQror:txx_xr;tyx_yzDxDyw_xw_yAx_aug offset_aug_xnc_augxnc_xAy_aug offset_aug_ync_augync_yry_Cr5_Ax_AyzhatRs&&&&&&&&&&&&& r1_solve_2d_fitpackrZsR ,,s C ,,s C%=  R%$!FL' 55D%=  R%$!FL' 55D Bw IIq"((BHHQK#.gsQx!|r3rg-q=)rr)rrinfmaxrrootr)r|rr}rr~r:rrr;rrrrfp0rrrrfpmsbracketftolrp_starfp_starC_starrs&&&&&&&&&&&&f&$$$ @r1_p_search_hit_srhsl bbars 'E R5DS1W~~.G q7{E "D!VWd%8A VVFFmG WWF 67 ""r3cb\VUu.uFqDRJNK up4'dWV\R4\R43#VwrVrxWV8dWx8g \R4h\P!W8W8*,4^,p \P!W8W8*,4^,p V P ^8XgV P ^8Xd \R4hW ,W,V\P !W4,\PV ,\PV ,3#uupi)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 < yez$bbox excludes all samples in x or y.)r!slicerrwhererix_s_) r)r*Zbboxbboxixbxeybyeixiys &&&& r1_apply_bbox_gridrs2 t ,teTMt ,--Qd U4[00NBB G@AA 17qw' ( +B 17qw' ( +B ww!|rww!|?@@ 5!%266">*BEE"IruuRy @@ -s D,c\P!V4p\P!V4p\P!WWW4wrp \P!WWh4wrpVP 4p\ WV 4\ WV4V3#r@)rrtr data_matrixrjr>)r)r*rrrr:r;rrr|offset_xrr}offset_yrr~s&&&&&&& r1_build_design_matricesrsv ,,q/C ,,q/C!--aR=B$!--aR=B$ A t , t ,  r3cVf)\W,^,^V,^,4pM6V^V^,,8d"\RV: R^V^,, R24h^V^,,pW,^,p\P!V.V^,,V.V^,,,4pWtWV3#)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) z`nest` too small: nest = z < 2*(k+1) = .)rrrr)rXrrrWnestnminnmaxrs&&&&& r1_initialise_knotsrsH |1519acAg& !QU) 9$-1Q3yPQRS S a!e9D 519D B41:ac *+A D r3c V\,p VPp Wr, p W8Xd^p MSW, pW8d\W,V, 4MV ^,p\V ^,\ W^,^44p \ V 4FBp\ WW4pVP^,p W8d\W4V 3u#W8gK?W:3u# W:3#)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). ) TOLrintrMrrLrrr)r)rWrrrrrr5fpold residualsnplusaccnrdeltanpl1js&&&&&&&&&&& r1 _add_knotsrsZ c'C A 6D  y ,1Ks5<%'(U1WE!GSax34 5\ Q1 ( GGAJ 9q$e+ + 98O'* 8Or3r:r;rrnestxnestyrc0\WW)4wrp p V PV^,8gV PV^,8dA\RV RV RV^, RV^, RV P RV P R2 4h\V ^,f V ^,MV ^,4p\V ^,f V R,MV ^,4p\V ^,f V ^,MV ^,4p\V ^,f V R,MV ^,4pRpVR8XdeVfVe \R 4h\ W4p\ W4p\ WVVVW44wppp\ VVVVVVWVW4 wppp \VVVV3W434#\V PWW7R 7wpppp\V PVVWHR 7wppppRpR p\V4\V4,p Rp!Rp"Rp#\V 4Fp \ WVVVW44wppp\ VVVVVVV VVV V 4 wppp$\V4V8Xd\V4V8XdTp!VV8dMVR 8Xd9\WVVVVVVV\P!V$^,^R 7V"R 7 wpp"RpM7\WVVVVVVV\P!V$^,^R 7V#R 7 wpp#R p\V4V8d\V4V8dMTpK \V4V8Xd"\V4V8Xd\XVVX3W434#^p\VV4wp%p&p'\VV4wp(p)p*\V%V&V'4p%\V(V)V*4p(\ WVVVW44wppp\!VV%VV(VVVWVWVV!VVR7wp p+p,\V,VVV+3W434#)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`. z/Not enough samples inside bbox for degrees (kx=z, ky=z ). Need at least k+1 per axis: (z, z). Got (z).Nrzs == 0 is interpolation only)rr*r)rrrr5rrrr))rrr)rrrrkrrrr<rrrLrrrzrr>r)-r)r*rr:r;rrrrrx_fity_fitZ_fitrrrrrrorrr|r}r~C0r5nminxnmaxxnminynmaxyr last_axismpmrnplusxnplusyrDrx offset_dxnc_dxDry offset_dync_dyC_smfp_sms-&&&$$$$$$$ r1_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 increasingz7x dimension of z must have same number of elements as xz7y dimension of z must have same number of elements as yzs should be s >= 0.0z"bbox shape should be (4,), found: )r:r;rrrrr)) rrrkr&r!r"rrrr)r)r*rrr:r;rrs&&&$$$$$r1_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 ++r3)rHrHT)NNr@)rcnumpyrscipy.interpolate._ndbspliner scipy.interpolate._fitpack_reprorrrrr scipy.sparserscipy._lib._utilr r2r<r>rrrrrrrrrrrrrr3r1rsTj2,,"*G0T3<))X&PJX??DF# F#F#%'F#R%AP 0f tnK=K=K=!K= K=K= $K= K=\4+T4+a4+A4+4+B4+r3