o iit@sddlZddlZddlmZddlmZddlmZddlZ ddgZ ddgZ d Z ej e d d Zd Zej ed d de DdZdZej ed dde DdZddZddZd$ddZddZddZd%ddZd%d d!ZGd"d#d#ZdS)&N)internal) get_typename) csr_matrixdoublezthrust::complexintz long longa$ #include extern "C" { __global__ void find_interval( const double* t, const double* x, long long* out, int k, int n, bool extrapolate, int total_x) { int idx = blockDim.x * blockIdx.x + threadIdx.x; if(idx >= total_x) { return; } double xp = *&x[idx]; double tb = *&t[k]; double te = *&t[n]; if(isnan(xp)) { out[idx] = -1; return; } if((xp < tb || xp > te) && !extrapolate) { out[idx] = -1; return; } int left = k; int right = n; int mid; bool found = false; while(left < right && !found) { mid = ((right + left) / 2); if(xp > *&t[mid]) { left = mid + 1; } else if (xp < *&t[mid]) { right = mid - 1; } else { found = true; } } int default_value = left - 1 < k ? k : left - 1; int result = found ? mid + 1 : default_value + 1; while(result != n && xp >= *&t[result]) { result++; } out[idx] = result - 1; } } )z -std=c++11)codeoptionsaK #include #include #define COMPUTE_LINEAR 0x1 template __global__ void d_boor( const double* t, const T* c, const int k, const int mu, const double* x, const long long* intervals, T* out, double* temp, int num_c, int mode, int num_x) { int idx = blockDim.x * blockIdx.x + threadIdx.x; if(idx >= num_x) { return; } double xp = *&x[idx]; long long interval = *&intervals[idx]; double* h = temp + idx * (2 * k + 1); double* hh = h + k + 1; int ind, j, n; double xa, xb, w; if(mode == COMPUTE_LINEAR && interval < 0) { for(j = 0; j < num_c; j++) { out[num_c * idx + j] = CUDART_NAN; } return; } /* * Perform k-m "standard" deBoor iterations * so that h contains the k+1 non-zero values of beta_{ell,k-m}(x) * needed to calculate the remaining derivatives. */ h[0] = 1.0; for (j = 1; j <= k - mu; j++) { for(int p = 0; p < j; p++) { hh[p] = h[p]; } h[0] = 0.0; for (n = 1; n <= j; n++) { ind = interval + n; xb = t[ind]; xa = t[ind - j]; if (xb == xa) { h[n] = 0.0; continue; } w = hh[n - 1]/(xb - xa); h[n - 1] += w*(xb - xp); h[n] = w*(xp - xa); } } /* * Now do m "derivative" recursions * to convert the values of beta into the mth derivative */ for (j = k - mu + 1; j <= k; j++) { for(int p = 0; p < j; p++) { hh[p] = h[p]; } h[0] = 0.0; for (n = 1; n <= j; n++) { ind = interval + n; xb = t[ind]; xa = t[ind - j]; if (xb == xa) { h[mu] = 0.0; continue; } w = ((double) j) * hh[n - 1]/(xb - xa); h[n - 1] -= w; h[n] = w; } } if(mode != COMPUTE_LINEAR) { return; } // Compute linear combinations for(j = 0; j < num_c; j++) { out[num_c * idx + j] = 0; for(n = 0; n < k + 1; n++) { out[num_c * idx + j] = ( out[num_c * idx + j] + c[(interval + n - k) * num_c + j] * ((T) h[n])); } } } cCg|]}d|dqS)zd_boor<>).0 type_namer r ]/home/ubuntu/veenaModal/venv/lib/python3.10/site-packages/cupyx/scipy/interpolate/_bspline.py sr)rrname_expressionsa= #include template __global__ void compute_design_matrix( const int k, const long long* intervals, double* bspline_basis, double* data, U* indices, int num_intervals) { int idx = blockDim.x * blockIdx.x + threadIdx.x; if(idx >= num_intervals) { return; } long long interval = *&intervals[idx]; double* work = bspline_basis + idx * (2 * k + 1); for(int j = 0; j <= k; j++) { int m = (k + 1) * idx + j; data[m] = work[j]; indices[m] = (U) (interval - k + j); } } cCr )zcompute_design_matrixdd|D}d|}|r|d|dn|}||}|S)NcSsg|]}t|jqSr )rdtype)r argr r rrsz$_get_module_func..z, Return np.complex128 for complex dtypes, np.float64 otherwise.)cupy issubdtypecomplexfloatingcomplex_float_rr r r _get_dtypesr%FcCs@t|}t|j}|j|dd}|rt|std|S)z~Convert the input into a C contiguous float array. NB: Upcasts half- and single-precision floats to double precision. F)copyz$Array must not contain infs or nans.)rascontiguousarrayr%rastypeisfiniteall ValueError)x check_finitedtypr r r_as_float_arrays  r/c Cs|jd|d}tj|tjd}ttd} | |jddddfd|||||||jdftt|jdd} t |jdd|d} tt d |} | |jddddfd|||||||| | d|jdf dS) a1 Evaluate a spline in the B-spline basis. Parameters ---------- t : ndarray, shape (n+k+1) knots c : ndarray, shape (n, m) B-spline coefficients xp : ndarray, shape (s,) Points to evaluate the spline at. nu : int Order of derivative to evaluate. extrapolate : int, optional Whether to extrapolate to ouf-of-bounds points, or to return NaNs. out : ndarray, shape (s, m) Computed values of the spline at each of the input points. This argument is modified in-place. rr$ find_intervalr2Nd_boor) shaper empty_likeint64rINTERVAL_MODULErnpprodempty D_BOOR_MODULE) tckxpnu extrapolateoutn intervalsinterval_kernelnum_ctemp d_boor_kernelr r r_evaluate_splines  rKc Cs(|jd|d}tj|tjd}ttd}||jddddfd|||||||jdft|jdd|d}ttd|} | |jddddfd|d |d||d |dd|jdf tj|jd|dtj d} tt d |} | |jddddfd|||| ||jdf| |fS) a Returns a design matrix in CSR format. Note that only indices is passed, but not indptr because indptr is already precomputed in the calling Python function design_matrix. Parameters ---------- x : array_like, shape (n,) Points to evaluate the spline at. t : array_like, shape (nt,) Sorted 1D array of knots. k : int B-spline degree. extrapolate : bool, optional Whether to extrapolate to ouf-of-bounds points. indices : ndarray, shape (n * (k + 1),) Preallocated indices of the final CSR array. Returns ------- data The data array of a CSR array of the b-spline design matrix. In each row all the basis elements are evaluated at the certain point (first row - x[0], ..., last row - x[-1]). indices The indices array of a CSR array of the b-spline design matrix. rr0r$r1r2r3r4r5Ncompute_design_matrix) r6rr7r8rr9r<r=zerosfloat64DESIGN_MAT_MODULE) r,r>r@rCindicesrErFrG bspline_basisrJdatadesign_mat_kernelr r r_make_design_matrixs.   rTr0c Cs(|dkr t|| S|\}}}||krtd||dftdfdt|jdd}zNt|D]G}||dd|d| d}||}|dd||dd|||}tj|t |f|jddf}|dd}|d8}q2Wnt y}ztd ||d}~ww|||fS) a Compute the spline representation of the derivative of a given spline Parameters ---------- tck : tuple of (t, c, k) Spline whose derivative to compute n : int, optional Order of derivative to evaluate. Default: 1 Returns ------- tck_der : tuple of (t2, c2, k2) Spline of order k2=k-n representing the derivative of the input spline. Notes ----- .. seealso:: :class:`scipy.interpolate.splder` See Also -------- splantider, splev, spalde rz@Order of derivative (n = %r) must be <= order of spline (k = %r)r4NNr0zIThe spline has internal repeated knots and is not differentiable %d times) splantiderr+slicelenr6rangerr_r:rMFloatingPointError) tckrEr>r?r@shjdter r rsplderFs4     "($   rccCs|dkr t|| S|\}}}tdfdt|jdd}t|D]T}||dd|d| d}||}tj|d| d|dd|d}tjtd|jdd||dg|df}tj|d||df}|d7}q#|||fS) a Compute the spline for the antiderivative (integral) of a given spline. Parameters ---------- tck : tuple of (t, c, k) Spline whose antiderivative to compute n : int, optional Order of antiderivative to evaluate. Default: 1 Returns ------- tck_ader : tuple of (t2, c2, k2) Spline of order k2=k+n representing the antiderivative of the input spline. See Also -------- splder, splev, spalde Notes ----- The `splder` function is the inverse operation of this function. Namely, ``splder(splantider(tck))`` is identical to `tck`, modulo rounding error. .. seealso:: :class:`scipy.interpolate.splantider` rNrUr0)axisr0rVr4) rcrYrZr6r[rcumsumr\rM)r^rEr>r?r@r_r`rar r rrXs    "(  rXc@seZdZdZdddZedddZedd Zedd d Z edd dZ d ddZ ddZ ddZ d!ddZd!ddZd"ddZdS)#BSplineaUnivariate spline in the B-spline basis. .. math:: S(x) = \sum_{j=0}^{n-1} c_j B_{j, k; t}(x) where :math:`B_{j, k; t}` are B-spline basis functions of degree `k` and knots `t`. Parameters ---------- t : ndarray, shape (n+k+1,) knots c : ndarray, shape (>=n, ...) spline coefficients k : int B-spline degree extrapolate : bool or 'periodic', optional whether to extrapolate beyond the base interval, ``t[k] .. t[n]``, or to return nans. If True, extrapolates the first and last polynomial pieces of b-spline functions active on the base interval. If 'periodic', periodic extrapolation is used. Default is True. axis : int, optional Interpolation axis. Default is zero. Attributes ---------- t : ndarray knot vector c : ndarray spline coefficients k : int spline degree extrapolate : bool If True, extrapolates the first and last polynomial pieces of b-spline functions active on the base interval. axis : int Interpolation axis. tck : tuple A read-only equivalent of ``(self.t, self.c, self.k)`` Notes ----- B-spline basis elements are defined via .. math:: B_{i, 0}(x) = 1, \textrm{if $t_i \le x < t_{i+1}$, otherwise $0$,} B_{i, k}(x) = \frac{x - t_i}{t_{i+k} - t_i} B_{i, k-1}(x) + \frac{t_{i+k+1} - x}{t_{i+k+1} - t_{i+1}} B_{i+1, k-1}(x) **Implementation details** - At least ``k+1`` coefficients are required for a spline of degree `k`, so that ``n >= k+1``. Additional coefficients, ``c[j]`` with ``j > n``, are ignored. - B-spline basis elements of degree `k` form a partition of unity on the *base interval*, ``t[k] <= x <= t[n]``. - Based on [1]_ and [2]_ .. seealso:: :class:`scipy.interpolate.BSpline` References ---------- .. [1] Tom Lyche and Knut Morken, Spline methods, http://www.uio.no/studier/emner/matnat/ifi/INF-MAT5340/v05/undervisningsmateriale/ .. [2] Carl de Boor, A practical guide to splines, Springer, 2001. TrcCs~t||_t||_tj|tjd|_|dkr||_ nt ||_ |jj d|jd}t ||jj}||_|dkrEt|j|d|_|dkrMtd|jjdkrWtd||jdkrjtdd|d|ft|jdkrxtd tt|j||ddkrtd t|jstd |jjdkrtd |jj d|krtd t|jj}tj|j|d|_dS)Nr$periodicrr0 Spline order cannot be negative.z$Knot vector must be one-dimensional.z$Need at least %d knots for degree %dr4z(Knots must be in a non-decreasing order.z!Need at least two internal knots.z#Knots should not have nans or infs.z,Coefficients must be at least 1-dimensional.z0Knots, coefficients and degree are inconsistent.)operatorindexr@rasarrayr?r'rNr>rCboolr6r_normalize_axis_indexndimrdmoveaxisr+diffanyrZuniquer)r*r%r)selfr>r?r@rCrdrErar r r__init__sB       zBSpline.__init__cCs0t|}||||_|_|_||_||_|S)zConstruct a spline without making checks. Accepts same parameters as the regular constructor. Input arrays `t` and `c` must of correct shape and dtype. )object__new__r>r?r@rCrd)clsr>r?r@rCrdrtr r rconstruct_fast/s zBSpline.construct_fastcCs|j|j|jfS)z@Equivalent to ``(self.t, self.c, self.k)`` (read-only). )r>r?r@rtr r rr^;sz BSpline.tckcCsbt|d}t|}tj|ddf|||ddf|f}t|}d||<|||||S)aoReturn a B-spline basis element ``B(x | t[0], ..., t[k+1])``. Parameters ---------- t : ndarray, shape (k+2,) internal knots extrapolate : bool or 'periodic', optional whether to extrapolate beyond the base interval, ``t[0] .. t[k+1]``, or to return nans. If 'periodic', periodic extrapolation is used. Default is True. Returns ------- basis_element : callable A callable representing a B-spline basis element for the knot vector `t`. Notes ----- The degree of the B-spline, `k`, is inferred from the length of `t` as ``len(t)-2``. The knot vector is constructed by appending and prepending ``k+1`` elements to internal knots `t`. .. seealso:: :class:`scipy.interpolate.BSpline` r4rr0rVg?)rZr/rr\ zeros_likery)rxr>rCr@r?r r r basis_elementAs , zBSpline.basis_elementFc Cst|d}t|d}|dkrt|}|dkrtd|jdks.t|dd|ddkr6td|d t|d |d krHtd |d |dkrh|j|d}|||||||||}d }n!|st|||kst |||j d|dkrtd |d |j d}||d}|t t j j krt j }nt j}t j||d|d}t jd|d|d|d|d} t|||||\} }t| || f|j d|j d|dfdS)a Returns a design matrix as a CSR format sparse array. Parameters ---------- x : array_like, shape (n,) Points to evaluate the spline at. t : array_like, shape (nt,) Sorted 1D array of knots. k : int B-spline degree. extrapolate : bool or 'periodic', optional Whether to extrapolate based on the first and last intervals or raise an error. If 'periodic', periodic extrapolation is used. Default is False. Returns ------- design_matrix : `csr_matrix` object Sparse matrix in CSR format where each row contains all the basis elements of the input row (first row = basis elements of x[0], ..., last row = basis elements x[-1]). Notes ----- In each row of the design matrix all the basis elements are evaluated at the certain point (first row - x[0], ..., last row - x[-1]). `nt` is a length of the vector of knots: as far as there are `nt - k - 1` basis elements, `nt` should be not less than `2 * k + 2` to have at least `k + 1` basis element. Out of bounds `x` raises a ValueError. .. note:: This method returns a `csr_matrix` instance as CuPy still does not have `csr_array`. .. seealso:: :class:`scipy.interpolate.BSpline` Trhrrir0NrVz2Expect t to be a 1-D sorted array_like, but got t=.r4zLength t is not enough for k=FzOut of bounds w/ x = r$)r6)r/rmr+ror:rrrZsizeminmaxr6riinfoint32r8r<arangerTr) rxr,r>r@rCrEnnz int_dtyperPindptrrRr r r design_matrixdsB ) ( $.  " zBSpline.design_matrixNc Cs<|dur|j}t|}|j|j}}tjt|tjd}|dkrF|jj |j d}|j|j ||j|j |j||j|j }d}tj t |t t|jjddf|jjd}||||||||jjdd}|jdkrtt|j}||||j|d||||jd}||}|S)a Evaluate a spline function. Parameters ---------- x : array_like points to evaluate the spline at. nu : int, optional derivative to evaluate (default is 0). extrapolate : bool or 'periodic', optional whether to extrapolate based on the first and last intervals or return nans. If 'periodic', periodic extrapolation is used. Default is `self.extrapolate`. Returns ------- y : array_like Shape is determined by replacing the interpolation axis in the coefficient array with the shape of `x`. Nr$rhr0Fr)rCrrlr6ror'ravelrNr>r~r@r<rZrr:r;r?r _evaluatereshaperdlistr[ transpose) rtr,rBrCx_shapex_ndimrErD dim_orderr r r__call__s4   &   zBSpline.__call__cCs4|jjjs |j|_|jjjs|j|_dSdSrU)r>flags c_contiguousr&r?rzr r r_ensure_c_contiguouss   zBSpline._ensure_c_contiguouscCs.t|j|j|jjdd|j||||dS)NrrV)rKr>r?rr6r@)rtrArBrCrDr r rrs zBSpline._evaluater0cCsn|j}t|jt|}|dkr"tj|t|f|jddf}t|j||jf|}|j ||j |j dS)al Return a B-spline representing the derivative. Parameters ---------- nu : int, optional Derivative order. Default is 1. Returns ------- b : BSpline object A new instance representing the derivative. See Also -------- splder, splantider rr0NrCrd) r?rZr>rr\rMr6rcr@ryrCrd)rtrBr?ctr^r r r derivatives$ zBSpline.derivativecCs|j}t|jt|}|dkr"tj|t|f|jddf}t|j||jf|}|j dkr4d}n|j }|j |||j dS)a Return a B-spline representing the antiderivative. Parameters ---------- nu : int, optional Antiderivative order. Default is 1. Returns ------- b : BSpline object A new instance representing the antiderivative. Notes ----- If antiderivative is computed and ``self.extrapolate='periodic'``, it will be set to False for the returned instance. This is done because the antiderivative is no longer periodic and its correct evaluation outside of the initially given x interval is difficult. See Also -------- splder, splantider rr0NrhFr) r?rZr>rr\rMr6rXr@rCryrd)rtrBr?rr^rCr r rantiderivatives$ zBSpline.antiderivativecCs4|dur|j}|d}||kr||}}d}|jj|jd}|dkr<|sr~r@ritemrrr<rr:r;r?r6rrZr\rMrXdivmodrlrNrKr)rtabrCsignrErDr?rtacakatsteperiodinterval n_periodsleftr,integralr r r integrateBsv  "$         zBSpline.integrate)Tr)TF)rNrerU)__name__ __module__ __qualname____doc__ru classmethodrypropertyr^r|rrrrrrrr r r rrgs" H.   "  [4  (rgrre)rjr cupy._corercupy._core._scalarrcupyx.scipy.sparsernumpyr:TYPES INT_TYPESINTERVAL_KERNEL RawModuler9 D_BOOR_KERNELr=DESIGN_MAT_KERNELrOrr%r/rKrTrcrXrgr r r rs@   6b   % 5 <6