o X۷i; @sdZddlmZddlZddlZddlZddlmZmZddl m Z m Z m Z m Z mZddZdZd Zd d Zd d ZddZddZddZdddddddddZdddddddddZddZddZGdddZd d!ZGd"d#d#Zd$d%d&Zd'd(Zdddddddd)d*d+Z dddddddd)d,d-Z!ddddddddd.d/d0Z"ddd1d1d2d2ddd3Z#Gd4d5d5Z$Gd6d7d7e$Z%Gd8d9d9e$Z&dS):a Replicate FITPACK's logic for smoothing spline functions and curves. Currently provides analogs of splrep and splprep python routines, i.e. curfit.f and parcur.f routines (the drivers are fpcurf.f and fppara.f, respectively) The Fortran sources are from https://github.com/scipy/scipy/blob/maintenance/1.11.x/scipy/interpolate/fitpack/ .. [1] P. Dierckx, "Algorithms for smoothing data with periodic and parametric splines, Computer Graphics and Image Processing", 20 (1982) 171-184. :doi:`10.1016/0146-664X(82)90043-0`. .. [2] P. Dierckx, "Curve and surface fitting with splines", Monographs on Numerical Analysis, Oxford University Press, 1993. .. [3] P. Dierckx, "An algorithm for smoothing, differentiation and integration of experimental data using spline functions", Journal of Computational and Applied Mathematics, vol. I, no 3, p. 165 (1975). https://doi.org/10.1016/0771-050X(75)90034-0 ) annotationsN)BSplinemake_interp_spline)fpback _not_a_knot _lsq_solve_qr QR_MODULE_get_module_funcc CsJt|}t|}|jdks|jdkrtd|d||jd}|jd}||d}|d|kr:|ksIntd|d|d|d|d |d|d|d krctd |d||d ||dd kr{td |d||d||||||dkrtd|d|d||ks|d |||dkrtd|d|dd|d|d}|d||dks|d |||d krt||jd}|d}||d}|d krd Std|dD]+} || } |d7}||} t|| k} | |dkrt||| | kr"t|qd S)z Check consistency of the data vector `x` and the knot vector `t`. Return None if inputs are consistent, raises a ValueError otherwise. z*Expect `x` and `t` be 1D sequences. Got x=z and t=rzNeed k+1 <= n-k-1 <= m. Got m=z, n=z and k=.Nz%First k knots must be ordered; got t=z$Last k knots must be ordered; got t=z'Internal knots must be distinct. Got t=zOut of bounds: x=z1Schoenberg-Whitney condition is violated with t =z and x =)cupyasarrayndim ValueErrorshapeanyrangeargmax) xtkmnnk1mesgllnk3jtjtlir$\/home/ubuntu/vllm_env/lib/python3.10/site-packages/cupyx/scipy/interpolate/_fitpack_repro.pyfpcheck#sP     $ ,(,  r&gMbP?c CsD|d}t|||||\}}}t|}t|||}t||||S)Nr )rrascontiguousarrayr_compute_residuals) ryrrww2_csplr$r$r%_get_residualss   r0cCs||djdd}||S)Nr r axissum)r,splxr*deltar$r$r%r)sr)cst|||| fddttdD}dd}tt|D]}||d}|||7<||d|8<q)|dd7<|fS)z)Split the knot interval into "runs". cs(g|]}||dqSr r3.0r#ix residualsr$r% s(z_split..r r r )r searchsortedrlen)rrrr<fpartscarriesr#carryr$r:r%_splits  rCc Cst||||\}}d}d}tt|D]}||d||dkr-|||kr-|}||}q|dkr6td||||ddd} || } t|| } tj|d| | || df} | S)apAdd a new knot. (Approximately) replicate FITPACK's logic: 1. split the `x` array into knot intervals, `t(j+k) <= x(i) <= t(j+k+1)` 2. find the interval with the maximum sum of residuals 3. insert a new knot into the middle of that interval. NB: a new knot is in fact an `x` value at the middle of the interval. So *the knots are a subset of `x`*. This routine is an analog of https://github.com/scipy/scipy/blob/v1.11.4/scipy/interpolate/fitpack/fpcurf.f#L190-L215 (cf _split function) and https://github.com/scipy/scipy/blob/v1.11.4/scipy/interpolate/fitpack/fpknot.f ig}Ô%Ir z4Internal error, please report it to CuPy developers.r N)rCrr?rrr>r_) rrrr<r@r;idx_max fpart_maxr# idx_newknotnew_knotidx_tt_newr$r$r%add_knots"$  rKcCstj|td}tj|td}|durtj|td}ntj|td}|jdkr/td|jd|dkr9td|jdksC|jdkrLtd |jd t|}|rd|jdkrctd |jd |d n|jdkrutd |jd |d |dddf}|jd|jdkrtd|jd|jd |jd|jdkrtd|jd|jd |jdks|dd|ddkrtdt |}|dkrtd||durt |}|durt |}|||||||fS)zACommon input validations for generate_knots and make_splrep. dtypeNr w.ndim=z not implemented yet.rzWeights must be non-negativer zy.ndim=z not supported (must be 1 or 2.)z$ != 2 not supported with parametric=r z$ != 1 not supported with parametric=z!Weights is incompatible: w.shape=z != zData is incompatible: x.shape=z and y.shape=r z(Expect `x` to be an ordered 1D sequence.z `s` must be non-negative. Got s=) rrfloat ones_likerrrboolroperatorindexminmax)rr*r+rsxbxe parametricr$r$r%_validate_inputssN    & rZrr+rWrXrrVnestc cs|dkr|dus |durtdt||}|VdSt|||||||t|dkd\}}}}}}}t||||||||dEdHdS)a Replicate FITPACK's constructing the knot vector. Parameters ---------- x, y : array_like The data points defining the curve ``y = f(x)``. w : array_like, optional Weights. xb : float, optional The boundary of the approximation interval. If None (default), is set to ``x[0]``. xe : float, optional The boundary of the approximation interval. If None (default), is set to ``x[-1]``. k : int, optional The spline degree. Default is cubic, ``k = 3``. s : float, optional The smoothing factor. Default is ``s = 0``. nest : int, optional Stop when at least this many knots are placed. Yields ------ t : ndarray Knot vectors with an increasing number of knots. The generator is finite: it stops when the smoothing critetion is satisfied, or when then number of knots exceeds the maximum value: the user-provided `nest` or `x.size + k + 1` --- which is the knot vector for the interpolating spline. Examples -------- Generate some noisy data and fit a sequence of LSQ splines: >>> import numpy as np >>> import matplotlib.pyplot as plt >>> from scipy.interpolate import make_lsq_spline, generate_knots >>> rng = np.random.default_rng(12345) >>> x = np.linspace(-3, 3, 50) >>> y = np.exp(-x**2) + 0.1 * rng.standard_normal(size=50) >>> knots = list(generate_knots(x, y, s=1e-10)) >>> for t in knots[::3]: ... spl = make_lsq_spline(x, y, t) ... xs = xs = np.linspace(-3, 3, 201) ... plt.plot(xs, spl(xs), '-', label=f'n = {len(t)}', lw=3, alpha=0.7) >>> plt.plot(x, y, 'o', label='data') >>> plt.plot(xs, np.exp(-xs**2), '--') >>> plt.legend() Note that increasing the number of knots make the result follow the data more and more closely. Also note that a step of the generator may add multiple knots: >>> [len(t) for t in knots] [8, 9, 10, 12, 16, 24, 40, 48, 52] Notes ----- The routine generates successive knots vectors of increasing length, starting from ``2*(k+1)`` to ``len(x) + k + 1``, trying to make knots more dense in the regions where the deviation of the LSQ spline from data is large. When the maximum number of knots, ``len(x) + k + 1`` is reached (this happens when ``s`` is small and ``nest`` is large), the generator stops, and the last output is the knots for the interpolation with the not-a-knot boundary condition. Knots are located at data sites, unless ``k`` is even and the number of knots is ``len(x) + k + 1``. In that case, the last output of the generator has internal knots at Greville sites, ``(x[1:] + x[:-1]) / 2``. rNzs == 0 is interpolation onlyr rYr[)rrrZrr_generate_knots_impl) rr*r+rWrXrrVr\rr$r$r%generate_knotssL r_c cs|t}|j} |durt| |dd|d}n|d|dkr1td|dd|ddd|d} | |d} tj|g|d|g|dtd} | jd} d }d }t| D]}| V|}t ||| ||d }| }||}t ||ks|dkrdS| | krd}n!||}||krt |||n|d}t |dt||dd}t|D]8}t|| ||} | jd} | | krt||} | VdS| |kr| VdS||dkrt ||| ||d }qq]dS) Nr r r`nest` too small: nest= < 2*(k+1) = r rLrr+)TOLsizerUrrrrOrrr0r4absintrTrKr)rr*r+rWrXrrVr\accrnminnmaxrrfpfpoldr-r<fpmsnplusr6npl1r r$r$r%r^osV  &       r^cCs>d}t|dD]}|||kr||||||9}q|S)Ng?r )r)rr#r rresrVr$r$r%prodds  rqc Cs|jd}|||d||}|d|d}tj|d|dftd}t|dD],}||d}t|dD]}||} || |d|| t|| |||||f<q:q,||||9}tddt|dD} ||d} || | fS)a{Discontinuity matrix. Matrix elements are jumps of k-th derivatives of b-splines at internal knots. See Eqs. (9)-(10) of Ref. [1], or, equivalently, Eq. (3.43) of Ref. [2]. This routine assumes internal knots are all simple (have multiplicity =1). Parameters ---------- t : ndarray, 1D, shape(n,) Knots. k : int The spline degree Returns ------- disc : ndarray, shape(n-2*k-1, k+2) The jumps of the k-th derivatives of b-splines at internal knots, ``t[k+1], ...., t[n-k-1]``. Notes ----- The normalization here follows FITPACK: (https://github.com/scipy/scipy/blob/maintenance/1.11.x/scipy/interpolate/fitpack/fpdisc.f#L36) The k-th derivative jumps are multiplied by a factor:: (delta / nrint)**k / k! where ``delta`` is the length of the interval spanned by internal knots, and ``nrint`` is one less the number of internal knots (i.e., the number of subintervals between them). References ---------- .. [1] Paul Dierckx, Algorithms for smoothing data with periodic and parametric splines, Computer Graphics and Image Processing, vol. 20, p. 171 (1982). :doi:`10.1016/0146-664X(82)90043-0` .. [2] Tom Lyche and Knut Morken, Spline methods, http://www.uio.no/studier/emner/matnat/ifi/INF-MAT5340/v05/undervisningsmateriale/ rr r rLcSsg|]}|qSr$r$r8r$r$r%r=#szdisc..)rremptyrOrrqarray) rrrr6nrintmatrjjr iir#offsetncr$r$r%discs / 0   rzc@s*eZdZdZddddddZddZdS) Fa The r.h.s. of ``f(p) = s``. Given scalar `p`, we solve the system of equations in the LSQ sense: | A | @ | c | = | y | | B / p | | 0 | | 0 | where `A` is the matrix of b-splines and `b` is the discontinuity matrix (the jumps of the k-th derivatives of b-spline basis elements at knots). Since we do that repeatedly while minimizing over `p`, we QR-factorize `A` only once and update the QR factorization only of the `B` rows of the augmented matrix |A, B/p|. The system of equations is Eq. (15) Ref. [1]_, the strategy and implementation follows that of FITPACK, see specific links below. References ---------- [1] P. Dierckx, Algorithms for Smoothing Data with Periodic and Parametric Splines, COMPUTER GRAPHICS AND IMAGE PROCESSING vol. 20, pp 171-184 (1982). https://doi.org/10.1016/0146-664X(82)90043-0 N)RYcCs||_||_||_||_|durtj|tdn|}|jdkr'td|jd||_ ||_ |jdkr;td|jdt ||\} } } |durV|durVt |||||\}}} |j d|d} |d}|j d|krztd |j dd |dd tj| j d|j dftd}tj|d| |f|_tj| | j d|jdftd}|d| ddf|d| d|f<||_tjtj| tjd| f|_| |_| |_dS) NrLr rNz != 1.r z$F: expected y.ndim == 2, got y.ndim=z instead.rzInternal error: R.shape[1]=z != k+1=r )rr*rrrrPrOrrr+rVrzrrzerosrDYYAAarangeintprxryb)selfrr*rrrVr+r|r}rb_offsetb_ncr-rynzzrr$r$r%__init__Cs8    "$ z F.__init__c Cs|j}|j}|j}|j|||dddf<|j}ttd}|dd||jd|jd||||jd|ft |||}t |j ||j }t |jd||j|j} | } ||_| |jS)N qr_reducer7rr r )rcopyrxryrrr rrrrrrr)r+rr*r4r/rV) rpABrxryQYrr.r/r<rkr$r$r%__call__us&       z F.__call__N)__name__ __module__ __qualname____doc__rrr$r$r$r%r{(s 2r{c Cs|||}|||}|||}|tjkr"|||| |S||||||||| ||||||S)aThe root of r(p) = (u*p + v) / (p + w) given three points and values, (p1, f2), (p2, f2) and (p3, f3). The FITPACK analog adjusts the bounds, and we do not https://github.com/scipy/scipy/blob/maintenance/1.11.x/scipy/interpolate/fitpack/fprati.f NB: FITPACK uses p < 0 to encode p=infinity. We just use the infinity itself. Since the bracket is ``p1 <= p2 <= p3``, ``p3`` can be infinite (in fact, this is what the minimizer starts with, ``p3=inf``). )rinf) p1f1p2f2p3f3h1h2h3r$r$r%fpratis   >rc@seZdZddZdS)BunchcKs|jjdi|dS)Nr$)__dict__update)rkwargsr$r$r%rszBunch.__init__N)rrrrr$r$r$r%rs ra9A theoretically impossible result was found during the iteration process for finding a smoothing spline with fp = s. probably causes : s too small. there is an approximation returned but the corresponding weighted sum of squared residuals does not satisfy the condition abs(fp-s)/s < tol. aHThe maximal number of iterations maxit (set to 20 by the program) allowed for finding a smoothing spline with fp=s has been reached. probably causes : s too small there is an approximation returned but the corresponding weighted sum of squared residuals does not satisfy the condition abs(fp-s)/s < tol. )r rcCs|d}d}d}d\}}|\\} } \} } |} ttD]}| || }}t||kr-d\}}nz|dkrR|| |krL|} |} | |} | | krK| |||} q|dkrRd}|dkr|| ||krv|} |} | |} | tjkru| | kru||| |} q|dkr|d}| |ks|| krd\}}nt| | ||| | } |dkr||} } q||} } qd \}}|dkrtjtt |d d t || ||d S) aSolve `f(p) = 0` using a rational function approximation. In a nutshell, since the function f(p) is known to be monotonically decreasing, we - maintain the bracket (p1, f1), (p2, f2) and (p3, f3) - at each iteration step, approximate f(p) by a rational function r(p) = (u*p + v) / (p + w) and make a step to p_new to the root of f(p): r(p_new) = 0. The coefficients u, v and w are found from the bracket values, p1..3 and f1...3 The algorithm and implementation follows https://github.com/scipy/scipy/blob/maintenance/1.11.x/scipy/interpolate/fitpack/fpcurf.f#L229 and https://github.com/scipy/scipy/blob/maintenance/1.11.x/scipy/interpolate/fitpack/fppara.f#L290 Note that the latter is for parametric splines and the former is for 1D spline functions. The minimization is identical though [modulo a summation over the dimensions in the computation of f(p)], so we reuse the minimizer for both d=1 and d>1. g?g?g{Gz?)rr)rTrr )r F)rFr ) stacklevel) convergedroot iterationsier) rMAXITrfrrrwarningswarnRuntimeWarning_iermesgr)fp0bracketrhcon1con9con4ich1ich3rrrrritrrrrr$r$r% root_ratisT      rr+rWrXrrVrr\c Cs|t} |j} |durt| |dd|d}n|d|dkr0td|dd|dd|dur8td|durOt||||||||d } t| d }nt||||jd d|dkrzt|||||\} } } d }t |d }t || ||fSt|||||\}}} |jd |d}||ddd f }t |||||d}| }||}t ||t |g|d|g|d||}| }||}d |ft j|ff}t||||||||d}t|||| }|j|fS)z3Shared infra for make_splrep and make_splprep. Nr r rr`rar zEither supply `t` or `nest`.)r+rrVrWrXr\r rrrc)rrVr+r|r})rdrerUrr^listr&rrrrr4r0rrsrr{rr/)rr*r+rWrXrrVrr\rhrgenr-r.rrpr|r}ryrr<rkfpinffp0rrr$r$r%_make_splrep_impl!sD  . rc Cs|dkr#|dus|dus|durtdtdd} | t_t|||dSt|||||||dd\}}}}}}}t|||||||||d \} } | t_| jdddf| _| S) aFind the B-spline representation of a 1D function. Given the set of data points ``(x[i], y[i])``, determine a smooth spline approximation of degree ``k`` on the interval ``xb <= x <= xe``. Parameters ---------- x, y : array_like, shape (m,) The data points defining a curve ``y = f(x)``. w : array_like, shape (m,), optional Strictly positive 1D array of weights, of the same length as `x` and `y`. The weights are used in computing the weighted least-squares spline fit. If the errors in the y values have standard-deviation given by the vector ``d``, then `w` should be ``1/d``. Default is ``np.ones(m)``. xb, xe : float, optional The interval to fit. If None, these default to ``x[0]`` and ``x[-1]``, respectively. k : int, optional The degree of the spline fit. It is recommended to use cubic splines, ``k=3``, which is the default. Even values of `k` should be avoided, especially with small `s` values. s : float, optional The smoothing condition. The amount of smoothness is determined by satisfying the conditions:: sum((w * (g(x) - y))**2 ) <= s where ``g(x)`` is the smoothed fit to ``(x, y)``. The user can use `s` to control the tradeoff between closeness to data and smoothness of fit. Larger `s` means more smoothing while smaller values of `s` indicate less smoothing. Recommended values of `s` depend on the weights, `w`. If the weights represent the inverse of the standard deviation of `y`, then a good `s` value should be found in the range ``(m-sqrt(2*m), m+sqrt(2*m))`` where ``m`` is the number of datapoints in `x`, `y`, and `w`. Default is ``s = 0.0``, i.e. interpolation. t : array_like, optional The spline knots. If None (default), the knots will be constructed automatically. There must be at least ``2*k + 2`` and at most ``m + k + 1`` knots. nest : int, optional The target length of the knot vector. Should be between ``2*(k + 1)`` (the minimum number of knots for a degree-``k`` spline), and ``m + k + 1`` (the number of knots of the interpolating spline). The actual number of knots returned by this routine may be slightly larger than `nest`. Default is None (no limit, add up to ``m + k + 1`` knots). Returns ------- spl : a `BSpline` instance For `s=0`, ``spl(x) == y``. For non-zero values of `s` the `spl` represents the smoothed approximation to `(x, y)`, generally with fewer knots. See Also -------- generate_knots : is used under the hood for generating the knots make_splprep : the analog of this routine for parametric curves make_interp_spline : construct an interpolating spline (``s = 0``) make_lsq_spline : construct the least-squares spline given the knot vector splrep : a FITPACK analog of this routine References ---------- .. [1] P. Dierckx, "Algorithms for smoothing data with periodic and parametric splines, Computer Graphics and Image Processing", 20 (1982) 171-184. .. [2] P. Dierckx, "Curve and surface fitting with splines", Monographs on Numerical Analysis, Oxford University Press, 1993. Notes ----- This routine constructs the smoothing spline function, :math:`g(x)`, to minimize the sum of jumps, :math:`D_j`, of the ``k``-th derivative at the internal knots (:math:`x_b < t_i < x_e`), where .. math:: D_i = g^{(k)}(t_i + 0) - g^{(k)}(t_i - 0) Specifically, the routine constructs the spline function :math:`g(x)` which minimizes .. math:: \sum_i | D_i |^2 \to \mathrm{min} provided that .. math:: \sum_{j=1}^m (w_j \times (g(x_j) - y_j))^2 \leqslant s , where :math:`s > 0` is the input parameter. In other words, we balance maximizing the smoothness (measured as the jumps of the derivative, the first criterion), and the deviation of :math:`g(x_j)` from the data :math:`y_j` (the second criterion). Note that the summation in the second criterion is over all data points, and in the first criterion it is over the internal spline knots (i.e. those with ``xb < t[i] < xe``). The spline knots are in general a subset of data, see `generate_knots` for details. Also note the difference of this routine to `make_lsq_spline`: the latter routine does not consider smoothness and simply solves a least-squares problem .. math:: \sum w_j \times (g(x_j) - y_j)^2 \to \mathrm{min} for a spline function :math:`g(x)` with a _fixed_ knot vector ``t``. rNs==0 is for interpolation onlyr r)rFr]r)rr make_splrep_resrrZrr.) rr*r+rWrXrrVrr\rpr/r$r$r%r_s{  r)r+uubuerrVrr\c Cstj|dd}|dur7|ddddf|ddddfd} t| jdd}tjd||df}|dkrV|dusG|dusG|durKtdt||j|dd|fSt |||||||d d \}}}}}}}t |||||||||d \} } | j j} t | j | | jdd} | |fS) ag Find a smoothed B-spline representation of a parametric N-D curve. Given a list of N 1D arrays, `x`, which represent a curve in N-dimensional space parametrized by `u`, find a smooth approximating spline curve ``g(u)``. Parameters ---------- x : array_like, shape (m, ndim) Sampled data points representing the curve in ``ndim`` dimensions. The typical use is a list of 1D arrays, each of length ``m``. w : array_like, shape(m,), optional Strictly positive 1D array of weights. The weights are used in computing the weighted least-squares spline fit. If the errors in the `x` values have standard deviation given by the vector d, then `w` should be 1/d. Default is ``cupy.ones(m)``. u : array_like, optional An array of parameter values for the curve in the parametric form. If not given, these values are calculated automatically, according to:: v[0] = 0 v[i] = v[i-1] + distance(x[i], x[i-1]) u[i] = v[i] / v[-1] ub, ue : float, optional The end-points of the parameters interval. Default to ``u[0]`` and ``u[-1]``. k : int, optional Degree of the spline. Cubic splines, ``k=3``, are recommended. Even values of `k` should be avoided especially with a small ``s`` value. Default is ``k=3`` s : float, optional A smoothing condition. The amount of smoothness is determined by satisfying the conditions:: sum((w * (g(u) - x))**2) <= s, where ``g(u)`` is the smoothed approximation to ``x``. The user can use `s` to control the trade-off between closeness and smoothness of fit. Larger ``s`` means more smoothing while smaller values of ``s`` indicate less smoothing. Recommended values of ``s`` depend on the weights, ``w``. If the weights represent the inverse of the standard deviation of ``x``, then a good ``s`` value should be found in the range ``(m - sqrt(2*m), m + sqrt(2*m))``, where ``m`` is the number of data points in ``x`` and ``w``. t : array_like, optional The spline knots. If None (default), the knots will be constructed automatically. There must be at least ``2*k + 2`` and at most ``m + k + 1`` knots. nest : int, optional The target length of the knot vector. Should be between ``2*(k + 1)`` (the minimum number of knots for a degree-``k`` spline), and ``m + k + 1`` (the number of knots of the interpolating spline). The actual number of knots returned by this routine may be slightly larger than `nest`. Default is None (no limit, add up to ``m + k + 1`` knots). Returns ------- spl : a `BSpline` instance For `s=0`, ``spl(u) == x``. For non-zero values of ``s``, `spl` represents the smoothed approximation to ``x``, generally with fewer knots. u : ndarray The values of the parameters See Also -------- generate_knots : is used under the hood for generating the knots make_splrep : the analog of this routine 1D functions make_interp_spline : construct an interpolating spline (``s = 0``) make_lsq_spline : construct the least-squares spline given the knot vector splprep : a FITPACK analog of this routine Notes ----- Given a set of :math:`m` data points in :math:`D` dimensions, :math:`\vec{x}_j`, with :math:`j=1, ..., m` and :math:`\vec{x}_j = (x_{j; 1}, ..., x_{j; D})`, this routine constructs the parametric spline curve :math:`g_a(u)` with :math:`a=1, ..., D`, to minimize the sum of jumps, :math:`D_{i; a}`, of the ``k``-th derivative at the internal knots (:math:`u_b < t_i < u_e`), where .. math:: D_{i; a} = g_a^{(k)}(t_i + 0) - g_a^{(k)}(t_i - 0) Specifically, the routine constructs the spline function :math:`g(u)` which minimizes .. math:: \sum_i \sum_{a=1}^D | D_{i; a} |^2 \to \mathrm{min} provided that .. math:: \sum_{j=1}^m \sum_{a=1}^D (w_j \times (g_a(u_j) - x_{j; a}))^2 \leqslant s where :math:`u_j` is the value of the parameter corresponding to the data point :math:`(x_{j; 1}, ..., x_{j; D})`, and :math:`s > 0` is the input parameter. In other words, we balance maximizing the smoothness (measured as the jumps of the derivative, the first criterion), and the deviation of :math:`g(u_j)` from the data :math:`x_j` (the second criterion). Note that the summation in the second criterion is over all data points, and in the first criterion it is over the internal spline knots (i.e. those with ``ub < t[i] < ue``). The spline knots are in general a subset of data, see `generate_knots` for details. References ---------- .. [1] P. Dierckx, "Algorithms for smoothing data with periodic and parametric splines, Computer Graphics and Image Processing", 20 (1982) 171-184. .. [2] P. Dierckx, "Curve and surface fitting with splines", Monographs on Numerical Analysis, Oxford University Press, 1993. r r1Nr r rr)rr2Tr]r)rstacksqrtr4cumsumrDrrTrZrr.rrr)rr+rrrrrVrr\dpr/rpccspl1r$r$r% make_splpreps$, rr r )r extrapolater r~r raiserconstc@seZdZdZddgddddfddZed!dd Zd d Zd d Zd"ddZ d#ddZ ddZ ddZ ddZ ddZddZd$ddZd$dd ZdS)%UnivariateSplinea] 1-D smoothing spline fit to a given set of data points. Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. `s` specifies the number of knots by specifying a smoothing condition. Parameters ---------- x : (N,) array_like 1-D array of independent input data. Must be increasing; must be strictly increasing if `s` is 0. y : (N,) array_like 1-D array of dependent input data, of the same length as `x`. w : (N,) array_like, optional Weights for spline fitting. Must be positive. If `w` is None, weights are all 1. Default is None. bbox : (2,) array_like, optional 2-sequence specifying the boundary of the approximation interval. If `bbox` is None, ``bbox=[x[0], x[-1]]``. Default is None. k : int, optional Degree of the smoothing spline. ``k = 3`` is a cubic spline. Default is 3. s : float or None, optional Positive smoothing factor used to choose the number of knots. Number of knots will be increased until the smoothing condition is satisfied:: sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s However, because of numerical issues, the actual condition is:: abs(sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) - s) < 0.001 * s If `s` is None, `s` will be set as `len(w)` for a smoothing spline that uses all data points. If 0, spline will interpolate through all data points. This is equivalent to `InterpolatedUnivariateSpline`. Default is None. The user can use the `s` to control the tradeoff between closeness and smoothness of fit. Larger `s` means more smoothing while smaller values of `s` indicate less smoothing. Recommended values of `s` depend on the weights, `w`. If the weights represent the inverse of the standard-deviation of `y`, then a good `s` value should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m is the number of datapoints in `x`, `y`, and `w`. This means ``s = len(w)`` should be a good value if ``1/w[i]`` is an estimate of the standard deviation of ``y[i]``. ext : int or str, optional Controls the extrapolation mode for elements not in the interval defined by the knot sequence. * if ext=0 or 'extrapolate', return the extrapolated value. * if ext=1 or 'zeros', return 0 * if ext=2 or 'raise', raise a ValueError * if ext=3 or 'const', return the boundary value. Default is 0. See Also -------- scipy.interpolate.UnivariateSpline Nr rrc Cs|durtdtj|td}tj|td}|\}} ||_|r!|n|d|_| r*| n|d|_||_||_||_ ||_ |durCt |}| || dSNz*weighted spline fitting is not implementedrLrr )NotImplementedErrorrrrOext_xb_xe_x_y_w_kr?set_smoothing_factor _reset_class) rrr*r+bboxrrVrrWrXr$r$r%rs$  zUnivariateSpline.__init__cCs,||}||_||_||_||_||_dSr)__new___splr _residualrr)clsr/residualrXrWrrr$r$r% _from_spls  zUnivariateSpline._from_splcCsZ|jj}|dkr dS|dkr|tdS|dkr |tdS|dkr+|tdSdS)Nrr rr )rr _set_classInterpolatedUnivariateSplineLSQUnivariateSpline)rrr$r$r%rszUnivariateSpline._reset_classcCs&||_|jtttfvr||_dS dSr) _spline_class __class__rrr)rrr$r$r%rs  zUnivariateSpline._set_classc Cs|j|j|j|jf\}}}}|j|j}}t|||||||d|_tj|_||_ |dur6t j |j dt d}|dur>|jj}t||dddf|||d|_|dS)z Continue spline computation with the given smoothing factor s and with the knots found at the last call. This routine modifies the spline in place. )rr+rWrXrVNrrLrc)rrrrrrrrr_sronesrrOrr0r4rr) rrVrrr*r+rrWrXr$r$r%rs$ z%UnivariateSpline.set_smoothing_factorcCs|||}|dur|j}nt|}|durtd||dkrf|j|j}}|dkr8d|||k||kB<|S|dkrPt||k||kBrNtd|d|S|d krf|||||k<|||||k<|S) a Evaluate spline (or its nu-th derivative) at positions x. Parameters ---------- x : ndarray A 1-D array of points at which to return the value of the smoothed spline or its derivatives. Note: `x` can be unordered but the evaluation is more efficient if `x` is (partially) ordered. nu : int The order of derivative of the spline to compute. ext : int Controls the value returned for elements of `x` not in the interval defined by the knot sequence. * if ext=0 or 'extrapolate', return the extrapolated value. * if ext=1 or 'zeros', return 0 * if ext=2 or 'raise', raise a ValueError * if ext=3 or 'const', return the boundary value. The default value is 0, passed from the initialization of UnivariateSpline. NzUnknown extrapolation mode %s.rr rbr zOut of bounds x=z with ext='raise'.r)rr _extrap_modesgetrrrr)rrnurresultrWrXr$r$r%r+s&   zUnivariateSpline.__call__cCs|jj}|jj|| S)z Return positions of interior knots of the spline. Internally, the knot vector contains ``2*k`` additional boundary knots. )rrr)rrr$r$r% get_knots[szUnivariateSpline.get_knotscCs|jjS)zReturn spline coefficients.)rr.rr$r$r% get_coeffscszUnivariateSpline.get_coeffscCs|jS)zReturn weighted sum of squared residuals of the spline approx. This is equivalent to:: sum((w[i] * (y[i]-spl(x[i])))**2, axis=0 )rrr$r$r% get_residualgszUnivariateSpline.get_residualcCsD||jkr ||jkp||jko||jk}|rtdS|j||S)al Return definite integral of the spline between two given points. Parameters ---------- a : float Lower limit of integration. b : float Upper limit of integration. Returns ------- integral : float The value of the definite integral of the spline between limits. rb)rrrrsr integrate)rarcondr$r$r%integralps  zUnivariateSpline.integralcs*fddtjjdD}tj|S)aReturn all derivatives of the spline at the point x. Parameters ---------- x : float The point to evaluate the derivatives at. Returns ------- der : ndarray, shape(k+1,) Derivatives of the orders 0 to k. csg|]}|qSr$)r)r9rrrr$r%r=sz0UnivariateSpline.derivatives..r )rrrrrD)rrlstr$rr% derivativess  zUnivariateSpline.derivativesr cCs:|j|}|jdkr dn|j}tj|||j|j|jdS)ak Construct a new spline representing the derivative of this spline. Parameters ---------- n : int, optional Order of derivative to evaluate. Default: 1 Returns ------- spline : UnivariateSpline Spline of order k2=k-n representing the derivative of this spline. rr rrrWrX)r derivativerrrrrr)rrr/rr$r$r%rs zUnivariateSpline.derivativecCs(|j|}tj||j|j|j|jdS)aw Construct a new spline representing the antiderivative of this spline. Parameters ---------- n : int, optional Order of antiderivative to evaluate. Default: 1 Returns ------- spline : UnivariateSpline Spline of order k2=k+n representing the antiderivative of this spline. r)rantiderivativerrrrrr)rrr/r$r$r%rs  zUnivariateSpline.antiderivative)rr)rNr7)rrrrr classmethodrrrrrrrrrrrrr$r$r$r%rs >  0  rcs0eZdZdZddgdddffdd ZZS)ra! 1-D interpolating spline for a given set of data points. Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. Spline function passes through all provided points. Equivalent to `UnivariateSpline` with `s` = 0. Parameters ---------- x : (N,) array_like Input dimension of data points -- must be strictly increasing y : (N,) array_like input dimension of data points w : (N,) array_like, optional Weights for spline fitting. Must be positive. If None (default), weights are all 1. bbox : (2,) array_like, optional 2-sequence specifying the boundary of the approximation interval. If None (default), ``bbox=[x[0], x[-1]]``. k : int, optional Degree of the smoothing spline. Default is ``k = 3``, a cubic spline. ext : int or str, optional Controls the extrapolation mode for elements not in the interval defined by the knot sequence. * if ext=0 or 'extrapolate', return the extrapolated value. * if ext=1 or 'zeros', return 0 * if ext=2 or 'raise', raise a ValueError * if ext=3 of 'const', return the boundary value. The default value is 0. See Also -------- scipy.interpolate.InterpolatedUnivariateSpline Nr rrc stj||d||||ddS)Nr)rVr+rrr)superr)rrr*r+rrrrr$r%rsz%InterpolatedUnivariateSpline.__init__)rrrrr __classcell__r$r$rr%rs$&rc@s(eZdZdZddgdddfddZdS)ra 1-D spline with explicit internal knots. Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. `t` specifies the internal knots of the spline Parameters ---------- x : (N,) array_like Input dimension of data points -- must be increasing y : (N,) array_like Input dimension of data points t : (M,) array_like interior knots of the spline. Must be in ascending order and:: bbox[0] < t[0] < ... < t[-1] < bbox[-1] w : (N,) array_like, optional weights for spline fitting. Must be positive. If None (default), weights are all 1. bbox : (2,) array_like, optional 2-sequence specifying the boundary of the approximation interval. If None (default), ``bbox = [x[0], x[-1]]``. k : int, optional Degree of the smoothing spline. Default is `k` = 3, a cubic spline. ext : int or str, optional Controls the extrapolation mode for elements not in the interval defined by the knot sequence. * if ext=0 or 'extrapolate', return the extrapolated value. * if ext=1 or 'zeros', return 0 * if ext=2 or 'raise', raise a ValueError * if ext=3 of 'const', return the boundary value. The default value is 0. Raises ------ ValueError If the interior knots do not satisfy the Schoenberg-Whitney conditions See Also -------- scipy.interpolate.LSQUnivariateSpline Nr rrc Cs|durtdtj|td}tj|td}|\}} ||_|r!|n|d|_| r*| n|d|_||_||_||_ ||_ t |} t |||| | ||dSr)rrrrOrrrrrrrr?r&rr) rrr*rr+rrrrWrXrVr$r$r%rs$   zLSQUnivariateSpline.__init__)rrrrrr$r$r$r%rs/r)'r __future__rrrRrcupyx.scipy.interpolaterr!cupyx.scipy.interpolate._bspline2rrrrr r&rdrr0r)rCrKrZr_r^rqrzr{rrrrrrrrrrrr$r$r$r%s` a(7 \ gJo ` >   /+