[Scipy-svn] r4687 - in trunk/scipy/interpolate: . src
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Sep 4 15:41:11 EDT 2008
Author: oliphant
Date: 2008-09-04 14:41:11 -0500 (Thu, 04 Sep 2008)
New Revision: 4687
Added:
trunk/scipy/interpolate/src/
trunk/scipy/interpolate/src/__fitpack.h
trunk/scipy/interpolate/src/_fitpackmodule.c
trunk/scipy/interpolate/src/fitpack.pyf
trunk/scipy/interpolate/src/multipack.h
Removed:
trunk/scipy/interpolate/__fitpack.h
trunk/scipy/interpolate/_fitpackmodule.c
trunk/scipy/interpolate/fitpack.pyf
trunk/scipy/interpolate/multipack.h
Modified:
trunk/scipy/interpolate/setup.py
Log:
Move interpolate sources to sub-directory.
Deleted: trunk/scipy/interpolate/__fitpack.h
===================================================================
--- trunk/scipy/interpolate/__fitpack.h 2008-09-04 16:47:55 UTC (rev 4686)
+++ trunk/scipy/interpolate/__fitpack.h 2008-09-04 19:41:11 UTC (rev 4687)
@@ -1,1128 +0,0 @@
-/*
- Python-C wrapper of FITPACK (by P. Dierckx) (in netlib known as dierckx)
- Author: Pearu Peterson <pearu at ioc.ee>
- June 1.-4., 1999
- June 7. 1999
- $Revision$
- $Date$
- */
-
-/* module_methods:
- {"_curfit", fitpack_curfit, METH_VARARGS, doc_curfit},
- {"_spl_", fitpack_spl_, METH_VARARGS, doc_spl_},
- {"_splint", fitpack_splint, METH_VARARGS, doc_splint},
- {"_sproot", fitpack_sproot, METH_VARARGS, doc_sproot},
- {"_spalde", fitpack_spalde, METH_VARARGS, doc_spalde},
- {"_parcur", fitpack_parcur, METH_VARARGS, doc_parcur},
- {"_surfit", fitpack_surfit, METH_VARARGS, doc_surfit},
- {"_bispev", fitpack_bispev, METH_VARARGS, doc_bispev},
- {"_insert", fitpack_insert, METH_VARARGS, doc_insert},
- */
-/* link libraries: (one item per line)
- ddierckx
- */
-/* python files: (to be imported to Multipack.py)
- fitpack.py
- */
-#if defined(NO_APPEND_FORTRAN)
-#define CURFIT curfit
-#define PERCUR percur
-#define SPALDE spalde
-#define SPLDER splder
-#define SPLEV splev
-#define SPLINT splint
-#define SPROOT sproot
-#define PARCUR parcur
-#define CLOCUR clocur
-#define SURFIT surfit
-#define BISPEV bispev
-#define PARDER parder
-#define INSERT insert
-#else
-#define CURFIT curfit_
-#define PERCUR percur_
-#define SPALDE spalde_
-#define SPLDER splder_
-#define SPLEV splev_
-#define SPLINT splint_
-#define SPROOT sproot_
-#define PARCUR parcur_
-#define CLOCUR clocur_
-#define SURFIT surfit_
-#define BISPEV bispev_
-#define PARDER parder_
-#define INSERT insert_
-#endif
-void CURFIT(int*,int*,double*,double*,double*,double*,double*,int*,double*,int*,int*,double*,double*,double*,double*,int*,int*,int*);
-void PERCUR(int*,int*,double*,double*,double*,int*,double*,int*,int*,double*,double*,double*,double*,int*,int*,int*);
-void SPALDE(double*,int*,double*,int*,double*,double*,int*);
-void SPLDER(double*,int*,double*,int*,int*,double*,double*,int*,double*,int*);
-void SPLEV(double*,int*,double*,int*,double*,double*,int*,int*);
-double SPLINT(double*,int*,double*,int*,double*,double*,double*);
-void SPROOT(double*,int*,double*,double*,int*,int*,int*);
-void PARCUR(int*,int*,int*,int*,double*,int*,double*,double*,double*,double*,int*,double*,int*,int*,double*,int*,double*,double*,double*,int*,int*,int*);
-void CLOCUR(int*,int*,int*,int*,double*,int*,double*,double*,int*,double*,int*,int*,double*,int*,double*,double*,double*,int*,int*,int*);
-void SURFIT(int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*,double*,int*,int*,int*,double*,int*,double*,int*,double*,double*,double*,double*,int*,double*,int*,int*,int*,int*);
-void BISPEV(double*,int*,double*,int*,double*,int*,int*,double*,int*,double*,int*,double*,double*,int*,int*,int*,int*);
-void PARDER(double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,int*,double*,int*,double*,double*,int*,int*,int*,int*);
-void INSERT(int*,double*,int*,double*,int*,double*,double*,int*,double*,int*,int*);
-
-/* Note that curev, cualde need no interface. */
-
-static char doc_bispev[] = " [z,ier] = _bispev(tx,ty,c,kx,ky,x,y,nux,nuy)";
-static PyObject *fitpack_bispev(PyObject *dummy, PyObject *args) {
- int nx,ny,kx,ky,mx,my,lwrk,*iwrk,kwrk,ier,lwa,mxy,nux,nuy;
- double *tx,*ty,*c,*x,*y,*z,*wrk,*wa = NULL;
- PyArrayObject *ap_x = NULL,*ap_y = NULL,*ap_z = NULL,*ap_tx = NULL,\
- *ap_ty = NULL,*ap_c = NULL;
- PyObject *x_py = NULL,*y_py = NULL,*c_py = NULL,*tx_py = NULL,*ty_py = NULL;
- if (!PyArg_ParseTuple(args, "OOOiiOOii",&tx_py,&ty_py,&c_py,&kx,&ky,
- &x_py,&y_py,&nux,&nuy))
- return NULL;
- ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1);
- ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, PyArray_DOUBLE, 0, 1);
- ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1);
- ap_tx = (PyArrayObject *)PyArray_ContiguousFromObject(tx_py, PyArray_DOUBLE, 0, 1);
- ap_ty = (PyArrayObject *)PyArray_ContiguousFromObject(ty_py, PyArray_DOUBLE, 0, 1);
- if (ap_x == NULL || ap_y == NULL || ap_c == NULL || ap_tx == NULL \
- || ap_ty == NULL) goto fail;
- x = (double *) ap_x->data;
- y = (double *) ap_y->data;
- c = (double *) ap_c->data;
- tx = (double *) ap_tx->data;
- ty = (double *) ap_ty->data;
- nx = ap_tx->dimensions[0];
- ny = ap_ty->dimensions[0];
- mx = ap_x->dimensions[0];
- my = ap_y->dimensions[0];
- mxy = mx*my;
- ap_z = (PyArrayObject *)PyArray_FromDims(1,&mxy,PyArray_DOUBLE);
- z = (double *) ap_z->data;
- if (nux || nuy)
- lwrk = mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1);
- else
- lwrk = mx*(kx+1)+my*(ky+1);
- kwrk = mx+my;
- lwa = lwrk+kwrk;
- if ((wa = (double *)malloc(lwa*sizeof(double)))==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- wrk = wa;
- iwrk = (int *)(wrk+lwrk);
- if (nux || nuy)
- PARDER(tx,&nx,ty,&ny,c,&kx,&ky,&nux,&nuy,x,&mx,y,&my,z,wrk,&lwrk,iwrk,&kwrk,&ier);
- else
- BISPEV(tx,&nx,ty,&ny,c,&kx,&ky,x,&mx,y,&my,z,wrk,&lwrk,iwrk,&kwrk,&ier);
-
- if (wa) free(wa);
- Py_DECREF(ap_x);
- Py_DECREF(ap_y);
- Py_DECREF(ap_c);
- Py_DECREF(ap_tx);
- Py_DECREF(ap_ty);
- return Py_BuildValue("Ni",PyArray_Return(ap_z),ier);
- fail:
- if (wa) free(wa);
- Py_XDECREF(ap_x);
- Py_XDECREF(ap_y);
- Py_XDECREF(ap_z);
- Py_XDECREF(ap_c);
- Py_XDECREF(ap_tx);
- Py_XDECREF(ap_ty);
- return NULL;
-}
-
-static char doc_surfit[] = " [tx,ty,c,o] = _surfit(x,y,z,w,xb,xe,yb,ye,kx,ky,iopt,s,eps,tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)";
-static PyObject *fitpack_surfit(PyObject *dummy, PyObject *args) {
- int iopt,m,kx,ky,nxest,nyest,nx,ny,lwrk1,lwrk2,*iwrk,kwrk,ier,lwa,nxo,nyo,\
- i,lc,lcest,nmax;
- double *x,*y,*z,*w,xb,xe,yb,ye,s,*tx,*ty,*c,fp,*wrk1,*wrk2,*wa = NULL,eps;
- PyArrayObject *ap_x = NULL,*ap_y = NULL,*ap_z,*ap_w = NULL,\
- *ap_tx = NULL,*ap_ty = NULL,*ap_c = NULL;
- PyArrayObject *ap_wrk = NULL;
- PyObject *x_py = NULL,*y_py = NULL,*z_py = NULL,*w_py = NULL,\
- *tx_py = NULL,*ty_py = NULL;
- PyObject *wrk_py=NULL;
- nx=ny=ier=nxo=nyo=0;
- if (!PyArg_ParseTuple(args, "OOOOddddiiiddOOiiOii",\
- &x_py,&y_py,&z_py,&w_py,&xb,&xe,\
- &yb,&ye,&kx,&ky,&iopt,&s,&eps,&tx_py,&ty_py,&nxest,&nyest,\
- &wrk_py,&lwrk1,&lwrk2)) return NULL;
- ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1);
- ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, PyArray_DOUBLE, 0, 1);
- ap_z = (PyArrayObject *)PyArray_ContiguousFromObject(z_py, PyArray_DOUBLE, 0, 1);
- ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, PyArray_DOUBLE, 0, 1);
- ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, PyArray_DOUBLE, 0, 1);
- /*ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, PyArray_INT, 0, 1);*/
- if (ap_x == NULL || ap_y == NULL || ap_z == NULL || ap_w == NULL \
- || ap_wrk == NULL) goto fail;
- x = (double *) ap_x->data;
- y = (double *) ap_y->data;
- z = (double *) ap_z->data;
- w = (double *) ap_w->data;
- m = ap_x->dimensions[0];
- nmax=nxest;
- if (nmax<nyest) nmax=nyest;
- lcest=(nxest-kx-1)*(nyest-ky-1);
- kwrk=m+(nxest-2*kx-1)*(nyest-2*ky-1);
- lwa = 2*nmax+lcest+lwrk1+lwrk2+kwrk;
- if ((wa = (double *)malloc(lwa*sizeof(double)))==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- tx = wa;
- ty = tx + nmax;
- c = ty + nmax;
- wrk1 = c + lcest;
- iwrk = (int *)(wrk1 + lwrk1);
- wrk2 = (double *)(iwrk+kwrk);
- if (iopt) {
- ap_tx=(PyArrayObject *)PyArray_ContiguousFromObject(tx_py, PyArray_DOUBLE, 0, 1);
- ap_ty=(PyArrayObject *)PyArray_ContiguousFromObject(ty_py, PyArray_DOUBLE, 0, 1);
- if (ap_tx == NULL || ap_ty == NULL) goto fail;
- nx = nxo = ap_tx->dimensions[0];
- ny = nyo = ap_ty->dimensions[0];
- memcpy(tx,ap_tx->data,nx*sizeof(double));
- memcpy(ty,ap_ty->data,ny*sizeof(double));
- }
- if (iopt==1) {
- lc = (nx-kx-1)*(ny-ky-1);
- memcpy(wrk1,ap_wrk->data,lc*sizeof(double));
- /*memcpy(iwrk,ap_iwrk->data,n*sizeof(int));*/
- }
- SURFIT(&iopt,&m,x,y,z,w,&xb,&xe,&yb,&ye,&kx,&ky,&s,&nxest,&nyest,&nmax,&eps,&nx,tx,&ny,ty,c,&fp,wrk1,&lwrk1,wrk2,&lwrk2,iwrk,&kwrk,&ier);
- i=0;
- while ((ier>10) && (i++<5)) {
- lwrk2=ier;
- if ((wrk2 = (double *)malloc(lwrk2*sizeof(double)))==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- SURFIT(&iopt,&m,x,y,z,w,&xb,&xe,&yb,&ye,&kx,&ky,&s,&nxest,&nyest,&nmax,&eps,&nx,tx,&ny,ty,c,&fp,wrk1,&lwrk1,wrk2,&lwrk2,iwrk,&kwrk,&ier);
- if (wrk2) free(wrk2);
- }
- if (ier==10) {
- PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
- goto fail;
- }
- lc = (nx-kx-1)*(ny-ky-1);
- Py_XDECREF(ap_tx);
- Py_XDECREF(ap_ty);
- ap_tx = (PyArrayObject *)PyArray_FromDims(1,&nx,PyArray_DOUBLE);
- ap_ty = (PyArrayObject *)PyArray_FromDims(1,&ny,PyArray_DOUBLE);
- ap_c = (PyArrayObject *)PyArray_FromDims(1,&lc,PyArray_DOUBLE);
- if (ap_tx == NULL || ap_ty == NULL || ap_c == NULL) goto fail;
- if ((iopt==0)||(nx>nxo)||(ny>nyo)) {
- Py_XDECREF(ap_wrk);
- ap_wrk = (PyArrayObject *)PyArray_FromDims(1,&lc,PyArray_DOUBLE);
- if (ap_wrk == NULL) goto fail;
- /*ap_iwrk = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_INT);*/
- }
- if(ap_wrk->dimensions[0]<lc) {
- Py_XDECREF(ap_wrk);
- ap_wrk = (PyArrayObject *)PyArray_FromDims(1,&lc,PyArray_DOUBLE);
- if (ap_wrk == NULL) goto fail;
- }
- memcpy(ap_tx->data,tx,nx*sizeof(double));
- memcpy(ap_ty->data,ty,ny*sizeof(double));
- memcpy(ap_c->data,c,lc*sizeof(double));
- memcpy(ap_wrk->data,wrk1,lc*sizeof(double));
- /*memcpy(ap_iwrk->data,iwrk,n*sizeof(int));*/
- if (wa) free(wa);
- Py_DECREF(ap_x);
- Py_DECREF(ap_y);
- Py_DECREF(ap_z);
- Py_DECREF(ap_w);
- return Py_BuildValue("NNN{s:N,s:i,s:d}",PyArray_Return(ap_tx),\
- PyArray_Return(ap_ty),PyArray_Return(ap_c),\
- "wrk",PyArray_Return(ap_wrk),\
- "ier",ier,"fp",fp);
- fail:
- if (wa) free(wa);
- Py_XDECREF(ap_x);
- Py_XDECREF(ap_y);
- Py_XDECREF(ap_z);
- Py_XDECREF(ap_w);
- Py_XDECREF(ap_tx);
- Py_XDECREF(ap_ty);
- Py_XDECREF(ap_wrk);
- /*Py_XDECREF(ap_iwrk);*/
- if (!PyErr_Occurred()) {
- PyErr_SetString(PyExc_ValueError, "An error occurred.");
- }
- return NULL;
-}
-
-
-static char doc_parcur[] = " [t,c,o] = _parcur(x,w,u,ub,ue,k,iopt,ipar,s,t,nest,wrk,iwrk,per)";
-static PyObject *fitpack_parcur(PyObject *dummy, PyObject *args) {
- int k,iopt,ipar,nest,*iwrk,idim,m,mx,n=0,no=0,nc,ier,lc,lwa,lwrk,i,per;
- double *x,*w,*u,*c,*t,*wrk,*wa=NULL,ub,ue,fp,s;
- PyObject *x_py = NULL,*u_py = NULL,*w_py = NULL,*t_py = NULL;
- PyObject *wrk_py=NULL,*iwrk_py=NULL;
- PyArrayObject *ap_x = NULL,*ap_u = NULL,*ap_w = NULL,*ap_t = NULL,*ap_c = NULL;
- PyArrayObject *ap_wrk = NULL,*ap_iwrk = NULL;
- if (!PyArg_ParseTuple(args, "OOOddiiidOiOOi",&x_py,&w_py,&u_py,&ub,&ue,\
- &k,&iopt,&ipar,&s,&t_py,&nest,&wrk_py,&iwrk_py,&per)) return NULL;
- ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1);
- ap_u = (PyArrayObject *)PyArray_ContiguousFromObject(u_py, PyArray_DOUBLE, 0, 1);
- ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, PyArray_DOUBLE, 0, 1);
- ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, PyArray_DOUBLE, 0, 1);
- ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, PyArray_INT, 0, 1);
- if (ap_x == NULL || ap_u == NULL || ap_w == NULL || ap_wrk == NULL || ap_iwrk == NULL) goto fail;
- x = (double *) ap_x->data;
- u = (double *) ap_u->data;
- w = (double *) ap_w->data;
- m = ap_w->dimensions[0];
- mx = ap_x->dimensions[0];
- idim = mx/m;
- if (per)
- lwrk=m*(k+1)+nest*(7+idim+5*k);
- else
- lwrk=m*(k+1)+nest*(6+idim+3*k);
- nc=idim*nest;
- lwa = nc+2*nest+lwrk;
- if ((wa = (double *)malloc(lwa*sizeof(double)))==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- t = wa;
- c = t + nest;
- wrk = c + nc;
- iwrk = (int *)(wrk + lwrk);
- if (iopt) {
- ap_t=(PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
- if (ap_t == NULL) goto fail;
- n = no = ap_t->dimensions[0];
- memcpy(t,ap_t->data,n*sizeof(double));
- }
- if (iopt==1) {
- memcpy(wrk,ap_wrk->data,n*sizeof(double));
- memcpy(iwrk,ap_iwrk->data,n*sizeof(int));
- }
- if (per)
- CLOCUR(&iopt,&ipar,&idim,&m,u,&mx,x,w,&k,&s,&nest,&n,t,&nc,\
- c,&fp,wrk,&lwrk,iwrk,&ier);
- else
- PARCUR(&iopt,&ipar,&idim,&m,u,&mx,x,w,&ub,&ue,&k,&s,&nest,&n,t,&nc,\
- c,&fp,wrk,&lwrk,iwrk,&ier);
- if (ier==10) goto fail;
- if (ier>0 && n==0) n=1;
- lc = (n-k-1)*idim;
- ap_t = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_DOUBLE);
- ap_c = (PyArrayObject *)PyArray_FromDims(1,&lc,PyArray_DOUBLE);
- if (ap_t == NULL || ap_c == NULL) goto fail;
- if ((iopt==0)||(n>no)) {
- ap_wrk = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_DOUBLE);
- ap_iwrk = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_INT);
- if (ap_wrk == NULL || ap_iwrk == NULL) goto fail;
- }
- memcpy(ap_t->data,t,n*sizeof(double));
- for (i=0;i<idim;i++)
- memcpy((double *) ap_c->data+i*(n-k-1),c+i*n,(n-k-1)*sizeof(double));
- memcpy(ap_wrk->data,wrk,n*sizeof(double));
- memcpy(ap_iwrk->data,iwrk,n*sizeof(int));
- if (wa) free(wa);
- Py_DECREF(ap_x);
- Py_DECREF(ap_w);
- return Py_BuildValue("NN{s:N,s:d,s:d,s:N,s:N,s:i,s:d}",PyArray_Return(ap_t),PyArray_Return(ap_c),"u",PyArray_Return(ap_u),"ub",ub,"ue",ue,"wrk",PyArray_Return(ap_wrk),"iwrk",PyArray_Return(ap_iwrk),"ier",ier,"fp",fp);
- fail:
- if (wa) free(wa);
- Py_XDECREF(ap_x);
- Py_XDECREF(ap_u);
- Py_XDECREF(ap_w);
- Py_XDECREF(ap_t);
- Py_XDECREF(ap_wrk);
- Py_XDECREF(ap_iwrk);
- return NULL;
-}
-
-static char doc_curfit[] = " [t,c,o] = _curfit(x,y,w,xb,xe,k,iopt,s,t,nest,wrk,iwrk,per)";
-static PyObject *fitpack_curfit(PyObject *dummy, PyObject *args) {
- int iopt,m,k,nest,n,lwrk,*iwrk,ier,lwa,lc,no=0,per;
- double *x,*y,*w,xb,xe,s,*t,*c,fp,*wrk,*wa = NULL;
- PyArrayObject *ap_x = NULL,*ap_y = NULL,*ap_w = NULL,*ap_t = NULL,*ap_c = NULL;
- PyArrayObject *ap_wrk = NULL,*ap_iwrk = NULL;
- PyObject *x_py = NULL,*y_py = NULL,*w_py = NULL,*t_py = NULL;
- PyObject *wrk_py=NULL,*iwrk_py=NULL;
- if (!PyArg_ParseTuple(args, "OOOddiidOiOOi",&x_py,&y_py,&w_py,&xb,&xe,\
- &k,&iopt,&s,&t_py,&nest,&wrk_py,&iwrk_py,&per)) return NULL;
- ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1);
- ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, PyArray_DOUBLE, 0, 1);
- ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, PyArray_DOUBLE, 0, 1);
- ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, PyArray_DOUBLE, 0, 1);
- ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, PyArray_INT, 0, 1);
- if (ap_x == NULL || ap_y == NULL || ap_w == NULL || ap_wrk == NULL || ap_iwrk == NULL) goto fail;
- x = (double *) ap_x->data;
- y = (double *) ap_y->data;
- w = (double *) ap_w->data;
- m = ap_x->dimensions[0];
- if (per) lwrk = m*(k+1) + nest*(8+5*k);
- else lwrk = m*(k+1) + nest*(7+3*k);
- lwa = 3*nest+lwrk;
- if ((wa = (double *)malloc(lwa*sizeof(double)))==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- t = wa;
- c = t + nest;
- wrk = c + nest;
- iwrk = (int *)(wrk + lwrk);
- if (iopt) {
- ap_t=(PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
- if (ap_t == NULL) goto fail;
- n = no = ap_t->dimensions[0];
- memcpy(t,ap_t->data,n*sizeof(double));
- }
- if (iopt==1) {
- memcpy(wrk,ap_wrk->data,n*sizeof(double));
- memcpy(iwrk,ap_iwrk->data,n*sizeof(int));
- }
- if (per)
- PERCUR(&iopt,&m,x,y,w,&k,&s,&nest,&n,t,c,&fp,wrk,&lwrk,iwrk,&ier);
- else
- CURFIT(&iopt,&m,x,y,w,&xb,&xe,&k,&s,&nest,&n,t,c,&fp,wrk,&lwrk,iwrk,&ier);
- if (ier==10) {
- PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
- goto fail;
- }
- lc = n-k-1;
- if (!iopt) {
- ap_t = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_DOUBLE);
- if (ap_t == NULL) goto fail;
- }
- ap_c = (PyArrayObject *)PyArray_FromDims(1,&lc,PyArray_DOUBLE);
- if (ap_c == NULL) goto fail;
- if ((iopt==0)||(n>no)) {
- Py_XDECREF(ap_wrk);
- Py_XDECREF(ap_iwrk);
- ap_wrk = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_DOUBLE);
- ap_iwrk = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_INT);
- if (ap_wrk == NULL || ap_iwrk == NULL) goto fail;
- }
- memcpy(ap_t->data,t,n*sizeof(double));
- memcpy(ap_c->data,c,lc*sizeof(double));
- memcpy(ap_wrk->data,wrk,n*sizeof(double));
- memcpy(ap_iwrk->data,iwrk,n*sizeof(int));
- if (wa) free(wa);
- Py_DECREF(ap_x);
- Py_DECREF(ap_y);
- Py_DECREF(ap_w);
- return Py_BuildValue("NN{s:N,s:N,s:i,s:d}",PyArray_Return(ap_t),PyArray_Return(ap_c),"wrk",PyArray_Return(ap_wrk),"iwrk",PyArray_Return(ap_iwrk),"ier",ier,"fp",fp);
- fail:
- if (wa) free(wa);
- Py_XDECREF(ap_x);
- Py_XDECREF(ap_y);
- Py_XDECREF(ap_w);
- Py_XDECREF(ap_t);
- Py_XDECREF(ap_wrk);
- Py_XDECREF(ap_iwrk);
- return NULL;
-}
-
-static char doc_spl_[] = " [y,ier] = _spl_(x,nu,t,c,k )";
-static PyObject *fitpack_spl_(PyObject *dummy, PyObject *args) {
- int n,nu,m,ier,k;
- double *x,*y,*t,*c,*wrk = NULL;
- PyArrayObject *ap_x = NULL,*ap_y = NULL,*ap_t = NULL,*ap_c = NULL;
- PyObject *x_py = NULL,*t_py = NULL,*c_py = NULL;
- if (!PyArg_ParseTuple(args, "OiOOi",&x_py,&nu,&t_py,&c_py,&k)) return NULL;
- ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1);
- ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
- ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1);
- if ((ap_x == NULL || ap_t == NULL || ap_c == NULL)) goto fail;
- x = (double *) ap_x->data;
- m = ap_x->dimensions[0];
- t = (double *) ap_t->data;
- c = (double *) ap_c->data;
- n = ap_t->dimensions[0];
- ap_y = (PyArrayObject *)PyArray_FromDims(1,&m,PyArray_DOUBLE);
- if (ap_y == NULL) goto fail;
- y = (double *) ap_y->data;
- if ((wrk = (double *)malloc(n*sizeof(double)))==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- if (nu)
- SPLDER(t,&n,c,&k,&nu,x,y,&m,wrk,&ier);
- else
- SPLEV(t,&n,c,&k,x,y,&m,&ier);
- if (wrk) free(wrk);
- Py_DECREF(ap_x);
- Py_DECREF(ap_c);
- Py_DECREF(ap_t);
- return Py_BuildValue("Ni",PyArray_Return(ap_y),ier);
- fail:
- if (wrk) free(wrk);
- Py_XDECREF(ap_x);
- Py_XDECREF(ap_c);
- Py_XDECREF(ap_t);
- return NULL;
-}
-
-static char doc_splint[] = " [aint,wrk] = _splint(t,c,k,a,b)";
-static PyObject *fitpack_splint(PyObject *dummy, PyObject *args) {
- int n,k;
- double *t,*c,*wrk = NULL,a,b,aint;
- PyArrayObject *ap_t = NULL,*ap_c = NULL;
- PyArrayObject *ap_wrk = NULL;
- PyObject *t_py = NULL,*c_py = NULL;
- if (!PyArg_ParseTuple(args, "OOidd",&t_py,&c_py,&k,&a,&b)) return NULL;
- ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
- ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1);
- if ((ap_t == NULL || ap_c == NULL)) goto fail;
- t = (double *) ap_t->data;
- c = (double *) ap_c->data;
- n = ap_t->dimensions[0];
- ap_wrk = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_DOUBLE);
- if (ap_wrk == NULL) goto fail;
- wrk = (double *) ap_wrk->data;
- aint = SPLINT(t,&n,c,&k,&a,&b,wrk);
- Py_DECREF(ap_c);
- Py_DECREF(ap_t);
- return Py_BuildValue("dN",aint,PyArray_Return(ap_wrk));
- fail:
- Py_XDECREF(ap_c);
- Py_XDECREF(ap_t);
- return NULL;
-}
-
-static char doc_sproot[] = " [z,ier] = _sproot(t,c,k,mest)";
-static PyObject *fitpack_sproot(PyObject *dummy, PyObject *args) {
- int n,k,mest,ier,m;
- double *t,*c,*z=NULL;
- PyArrayObject *ap_t = NULL,*ap_c = NULL;
- PyArrayObject *ap_z = NULL;
- PyObject *t_py = NULL,*c_py = NULL;
- if (!PyArg_ParseTuple(args, "OOii",&t_py,&c_py,&k,&mest)) return NULL;
- ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
- ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1);
- if ((ap_t == NULL || ap_c == NULL)) goto fail;
- t = (double *) ap_t->data;
- c = (double *) ap_c->data;
- n = ap_t->dimensions[0];
- if ((z = (double *)malloc(mest*sizeof(double)))==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- SPROOT(t,&n,c,z,&mest,&m,&ier);
- if (ier==10) m=0;
- ap_z = (PyArrayObject *)PyArray_FromDims(1,&m,PyArray_DOUBLE);
- if (ap_z == NULL) goto fail;
- memcpy(ap_z->data,z,m*sizeof(double));
- if (z) free(z);
- Py_DECREF(ap_c);
- Py_DECREF(ap_t);
- return Py_BuildValue("Ni",PyArray_Return(ap_z),ier);
- fail:
- if (z) free(z);
- Py_XDECREF(ap_c);
- Py_XDECREF(ap_t);
- return NULL;
-}
-
-static char doc_spalde[] = " [d,ier] = _spalde(t,c,k,x)";
-static PyObject *fitpack_spalde(PyObject *dummy, PyObject *args) {
- int n,k,k1,ier;
- double *t,*c,*d=NULL,x;
- PyArrayObject *ap_t = NULL,*ap_c = NULL,*ap_d = NULL;
- PyObject *t_py = NULL,*c_py = NULL;
- if (!PyArg_ParseTuple(args, "OOid",&t_py,&c_py,&k,&x)) return NULL;
- ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
- ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1);
- if ((ap_t == NULL || ap_c == NULL)) goto fail;
- t = (double *) ap_t->data;
- c = (double *) ap_c->data;
- n = ap_t->dimensions[0];
- k1=k+1;
- ap_d = (PyArrayObject *)PyArray_FromDims(1,&k1,PyArray_DOUBLE);
- if (ap_d == NULL) goto fail;
- d = (double *) ap_d->data;
- SPALDE(t,&n,c,&k1,&x,d,&ier);
- Py_DECREF(ap_c);
- Py_DECREF(ap_t);
- return Py_BuildValue("Ni",PyArray_Return(ap_d),ier);
- fail:
- Py_XDECREF(ap_c);
- Py_XDECREF(ap_t);
- return NULL;
-}
-
-static char doc_insert[] = " [tt,cc,ier] = _insert(iopt,t,c,k,x,m)";
-static PyObject *fitpack_insert(PyObject *dummy, PyObject*args) {
- int iopt, n, nn, k, nest, ier, m;
- double x;
- double *t, *c, *tt, *cc;
- PyArrayObject *ap_t = NULL, *ap_c = NULL, *ap_tt = NULL, *ap_cc = NULL;
- PyObject *t_py = NULL, *c_py = NULL;
- PyObject *ret = NULL;
- if (!PyArg_ParseTuple(args, "iOOidi",&iopt,&t_py,&c_py,&k, &x, &m)) return NULL;
- ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
- ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1);
- if (ap_t == NULL || ap_c == NULL) goto fail;
- t = (double *) ap_t->data;
- c = (double *) ap_c->data;
- n = ap_t->dimensions[0];
- nest = n + m;
- ap_tt = (PyArrayObject *)PyArray_FromDims(1,&nest,PyArray_DOUBLE);
- ap_cc = (PyArrayObject *)PyArray_FromDims(1,&nest,PyArray_DOUBLE);
- if (ap_tt == NULL || ap_cc == NULL) goto fail;
- tt = (double *) ap_tt->data;
- cc = (double *) ap_cc->data;
- for ( ; n < nest; n++) {
- INSERT(&iopt, t, &n, c, &k, &x, tt, &nn, cc, &nest, &ier);
- if (ier) break;
- t = tt;
- c = cc;
- }
- Py_DECREF(ap_c);
- Py_DECREF(ap_t);
- ret = Py_BuildValue("NNi",PyArray_Return(ap_tt),PyArray_Return(ap_cc),ier);
- return ret;
-
- fail:
- Py_XDECREF(ap_c);
- Py_XDECREF(ap_t);
- return NULL;
- }
-
-
-static void
-_deBoor_D(double *t, double x, int k, int ell, int m, double *result) {
- /* On completion the result array stores
- the k+1 non-zero values of beta^(m)_i,k(x): for i=ell, ell-1, ell-2, ell-k.
- Where t[ell] <= x < t[ell+1].
- */
- /* Implements a recursive algorithm similar to the original algorithm of
- deBoor.
- */
-
- double *hh = result + k + 1;
- double *h = result;
- double xb, xa, w;
- int ind, j, n;
-
- /* 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. */
-
- result[0] = 1.0;
- for (j=1; j<=k-m; j++) {
- memcpy(hh, h, j*sizeof(double));
- h[0] = 0.0;
- for (n=1; n<=j; n++) {
- ind = ell + 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-x);
- h[n] = w*(x-xa);
- }
- }
-
- /* Now do m "derivative" recursions */
- /* to convert the values of beta into the mth derivative */
- for (j=k-m+1; j<=k; j++) {
- memcpy(hh, h, j*sizeof(double));
- h[0] = 0.0;
- for (n=1; n<=j; n++) {
- ind = ell + n;
- xb = t[ind];
- xa = t[ind-j];
- if (xb == xa) {
- h[m] = 0.0;
- continue;
- }
- w = j*hh[n-1]/(xb-xa);
- h[n-1] -= w;
- h[n] = w;
- }
- }
-}
-
-
-/* Given a set of (N+1) samples: A default set of knots is constructed
- using the samples xk plus 2*(K-1) additional knots where
- K = max(order,1) and the knots are chosen so that distances
- are symmetric around the first and last samples: x_0 and x_N.
-
- There should be a vector of N+K coefficients for the spline
- curve in coef. These coefficients form the curve as
-
- s(x) = sum(c_j B_{j,K}(x), j=-K..N-1)
-
- The spline function is evaluated at all points xx.
- The approximation interval is from xk[0] to xk[-1]
- Any xx outside that interval is set automatically to 0.0
- */
-static char doc_bspleval[] = "y = _bspleval(xx,xk,coef,k,{deriv (0)})\n"
- "\n"
- "The spline is defined by the approximation interval xk[0] to xk[-1],\n"
- "the length of xk (N+1), the order of the spline, k, and \n"
- "the number of coeficients N+k. The coefficients range from xk_{-K}\n"
- "to xk_{N-1} inclusive and are all the coefficients needed to define\n"
- "an arbitrary spline of order k, on the given approximation interval\n"
- "\n"
- "Extra knot points are internally added using knot-point symmetry \n"
- "around xk[0] and xk[-1]";
-
-static PyObject *_bspleval(PyObject *dummy, PyObject *args) {
- int k,kk,N,i,ell,dk,deriv=0;
- PyObject *xx_py=NULL, *coef_py=NULL, *x_i_py=NULL;
- PyArrayObject *xx=NULL, *coef=NULL, *x_i=NULL, *yy=NULL;
- PyArrayIterObject *xx_iter;
- double *t=NULL, *h=NULL, *ptr;
- double x0, xN, xN1, arg, sp, cval;
- if (!PyArg_ParseTuple(args, "OOOi|i", &xx_py, &x_i_py, &coef_py, &k, &deriv))
- return NULL;
- if (k < 0) {
- PyErr_Format(PyExc_ValueError, "order (%d) must be >=0", k);
- return NULL;
- }
- if (deriv > k) {
- PyErr_Format(PyExc_ValueError, "derivative (%d) must be <= order (%d)",
- deriv, k);
- return NULL;
- }
- kk = k;
- if (k==0) kk = 1;
- dk = (k == 0 ? 0 : 1);
- x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
- coef = (PyArrayObject *)PyArray_FROMANY(coef_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
- xx = (PyArrayObject *)PyArray_FROMANY(xx_py, NPY_DOUBLE, 0, 0, NPY_ALIGNED);
- if (x_i == NULL || coef == NULL || xx == NULL) goto fail;
-
- N = PyArray_DIM(x_i,0)-1;
-
- if (PyArray_DIM(coef,0) < (N+k)) {
- PyErr_Format(PyExc_ValueError, "too few coefficients (have %d need at least %d)",
- PyArray_DIM(coef,0), N+k);
- goto fail;
- }
-
- /* create output values */
- yy = (PyArrayObject *)PyArray_EMPTY(xx->nd, xx->dimensions, NPY_DOUBLE, 0);
- if (yy == NULL) goto fail;
- /* create dummy knot array with new knots inserted at the end
- selected as mirror symmetric versions of the old knots
- */
- t = (double *)malloc(sizeof(double)*(N+2*kk-1));
- if (t==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- x0 = *((double *)PyArray_DATA(x_i));
- xN = *((double *)PyArray_DATA(x_i) + N);
- for (i=0; i<kk-1; i++) { /* fill in ends if kk > 1*/
- t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i,kk-1-i)));
- t[kk+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i,N-1-i)));
- }
- ptr = t + (kk-1);
- for (i=0; i<=N; i++) {
- *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)));
- }
-
- /* Create work array to hold computed non-zero values for
- the spline for a value of x.
- */
- h = (double *)malloc(sizeof(double)*(2*kk+1));
- if (h==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
-
- /* Determine the spline for each value of x */
- xx_iter = (PyArrayIterObject *)PyArray_IterNew((PyObject *)xx);
- if (xx_iter == NULL) goto fail;
- ptr = PyArray_DATA(yy);
-
- while(PyArray_ITER_NOTDONE(xx_iter)) {
- arg = *((double *)PyArray_ITER_DATA(xx_iter));
- if ((arg < x0) || (arg > xN)) {
- /* If we are outside the interpolation region,
- fill with zeros
- */
- *ptr++ = 0.0;
- }
- else {
- /* Find the interval that arg lies between in the set of knots
- t[ell] <= arg < t[ell+1] (last-knot use the previous interval) */
- xN1 = *((double *)PyArray_DATA(x_i) + N-1);
- if (arg >= xN1) {
- ell = N + kk - 2;
- }
- else {
- ell = kk-1;
- while ((arg > t[ell])) ell++;
- if (arg != t[ell]) ell--;
- }
-
- _deBoor_D(t, arg, k, ell, deriv, h);
-
- sp = 0.0;
- for (i=0; i<=k; i++) {
- cval = *((double *)(PyArray_GETPTR1(coef, ell-i+dk)));
- sp += cval*h[k-i];
- }
- *ptr++ = sp;
- }
- PyArray_ITER_NEXT(xx_iter);
- }
- Py_DECREF(xx_iter);
- Py_DECREF(x_i);
- Py_DECREF(coef);
- Py_DECREF(xx);
- free(t);
- free(h);
- return PyArray_Return(yy);
-
- fail:
- Py_XDECREF(xx);
- Py_XDECREF(coef);
- Py_XDECREF(x_i);
- Py_XDECREF(yy);
- if (t != NULL) free(t);
- if (h != NULL) free(h);
- return NULL;
-}
-
-
-/* Given a set of (N+1) sample positions:
- Construct the diagonals of the (N+1) x (N+K) matrix that is needed to find
- the coefficients of a spline fit of order K.
- Note that K>=2 because for K=0,1, the coefficients are just the
- sample values themselves.
-
- The equation that expresses the constraints is
-
- s(x_i) = sum(c_j B_{j,K}(x_i), j=-K..N-1) = w_i for i=0..N
-
- This is equivalent to
-
- w = B*c where c.T = [c_{-K}, c{-K+1}, ..., c_{N-1}] and
- w.T = [w_{0}, w_{1}, ..., w_{N}]
-
- Therefore B is an (N+1) times (N+K) matrix with entries
-
- B_{j,K}(x_i) for column j=-K..N-1
- and row i=0..N
-
- This routine takes the N+1 sample positions and the order k and
- constructs the banded constraint matrix B (with k non-zero diagonals)
-
- The returned array is (N+1) times (N+K) ready to be either used
- to compute a minimally Kth-order derivative discontinuous spline
- or to be expanded with an additional K-1 constraints to be used in
- an exact spline specification.
- */
-static char doc_bsplmat[] = "B = _bsplmat(order,xk)\n"
-"Construct the constraint matrix for spline fitting of order k\n"
-"given sample positions in xk.\n"
-"\n"
-"If xk is an integer (N+1), then the result is equivalent to\n"
-"xk=arange(N+1)+x0 for any value of x0. This produces the\n"
-"integer-spaced, or cardinal spline matrix a bit faster.";
-static PyObject *_bsplmat(PyObject *dummy, PyObject *args) {
- int k,N,i,numbytes,j, equal;
- npy_intp dims[2];
- PyObject *x_i_py=NULL;
- PyArrayObject *x_i=NULL, *BB=NULL;
- double *t=NULL, *h=NULL, *ptr;
- double x0, xN, arg;
- if (!PyArg_ParseTuple(args, "iO", &k, &x_i_py))
- return NULL;
- if (k < 2) {
- PyErr_Format(PyExc_ValueError, "order (%d) must be >=2", k);
- return NULL;
- }
-
- equal = 0;
- N = PySequence_Length(x_i_py);
- if (N == -1 && PyErr_Occurred()) {
- PyErr_Clear();
- N = PyInt_AsLong(x_i_py);
- if (N==-1 && PyErr_Occurred()) goto fail;
- equal = 1;
- }
- N -= 1;
-
- /* create output matrix */
- dims[0] = N+1;
- dims[1] = N+k;
- BB = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);
- if (BB == NULL) goto fail;
-
- t = (double *)malloc(sizeof(double)*(N+2*k-1));
- if (t==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
-
- /* Create work array to hold computed non-zero values for
- the spline for a value of x.
- */
- h = (double *)malloc(sizeof(double)*(2*k+1));
- if (h==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
-
- numbytes = k*sizeof(double);
-
- if (equal) { /* points equally spaced by 1 */
- /* we run deBoor's algorithm one time with artificially created knots
- Then, we keep copying the result to every row */
-
- /* Create knots at equally-spaced locations from -(K-1) to N+K-1 */
- ptr = t;
- for (i=-k+1; i<N+k; i++) *ptr++ = i;
- j = k-1;
- _deBoor_D(t, 0, k, k-1, 0, h);
- ptr = PyArray_DATA(BB);
- N = N+1;
- for (i=0; i<N; i++) {
- memcpy(ptr, h, numbytes);
- ptr += (N+k);
- }
- goto finish;
- }
-
- /* Not-equally spaced */
- x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
- if (x_i == NULL) return NULL;
-
- /* create dummy knot array with new knots inserted at the end
- selected as mirror symmetric versions of the old knots
- */
- x0 = *((double *)PyArray_DATA(x_i));
- xN = *((double *)PyArray_DATA(x_i) + N);
- for (i=0; i<k-1; i++) { /* fill in ends if k > 1*/
- t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i,k-1-i)));
- t[k+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i,N-1-i)));
- }
- ptr = t + (k-1);
- for (i=0; i<=N; i++) {
- *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)));
- }
-
-
- /* Determine the K+1 non-zero values of the spline and place them in the
- correct location in the matrix for each row (along the diagonals).
- In fact, the last member is always zero so only K non-zero values
- are present.
- */
- ptr = PyArray_DATA(BB);
- for (i=0,j=k-1; i<N; i++,j++) {
- arg = *((double *)PyArray_DATA(x_i) + i);
- _deBoor_D(t, arg, k, j, 0, h);
- memcpy(ptr, h, numbytes);
- ptr += (N+k+1); /* advance to next row shifted over one */
- }
- /* Last row is different the first coefficient is zero.*/
- _deBoor_D(t, xN, k, j-1, 0, h);
- memcpy(ptr, h+1, numbytes);
-
- finish:
- Py_XDECREF(x_i);
- free(t);
- free(h);
- return (PyObject *)BB;
-
- fail:
- Py_XDECREF(x_i);
- Py_XDECREF(BB);
- if (t != NULL) free(t);
- if (h != NULL) free(h);
- return NULL;
-}
-
-
-
-/* Given a set of (N+1) sample positions:
- Construct the (N-1) x (N+K) error matrix J_{ij} such that
-
- for i=1..N-1,
-
- e_i = sum(J_{ij}c_{j},j=-K..N-1)
-
- is the discontinuity of the Kth derivative at the point i in the spline.
-
- This routine takes the N+1 sample positions and the order k and
- constructs the banded matrix J
-
- The returned array is (N+1) times (N+K) ready to be either used
- to compute a minimally Kth-order derivative discontinuous spline
- or to be expanded with an additional K-1 constraints to be used in
- an exact reconstruction approach.
- */
-static char doc_bspldismat[] = "B = _bspldismat(order,xk)\n"
-"Construct the kth derivative discontinuity jump constraint matrix \n"
-"for spline fitting of order k given sample positions in xk.\n"
-"\n"
-"If xk is an integer (N+1), then the result is equivalent to\n"
-"xk=arange(N+1)+x0 for any value of x0. This produces the\n"
-"integer-spaced matrix a bit faster. If xk is a 2-tuple (N+1,dx)\n"
-"then it produces the result as if the sample distance were dx";
-static PyObject *_bspldismat(PyObject *dummy, PyObject *args) {
- int k,N,i,j, equal, m;
- npy_intp dims[2];
- PyObject *x_i_py=NULL;
- PyArrayObject *x_i=NULL, *BB=NULL;
- double *t=NULL, *h=NULL, *ptr, *dptr;
- double x0, xN, dx;
- if (!PyArg_ParseTuple(args, "iO", &k, &x_i_py))
- return NULL;
- if (k < 2) {
- PyErr_Format(PyExc_ValueError, "order (%d) must be >=2", k);
- return NULL;
- }
-
- equal = 0;
- N = PySequence_Length(x_i_py);
- if (N==2 || (N == -1 && PyErr_Occurred())) {
- PyErr_Clear();
- if (PyTuple_Check(x_i_py)) {
- /* x_i_py = (N+1, dx) */
- N = PyInt_AsLong(PyTuple_GET_ITEM(x_i_py, 0));
- dx = PyFloat_AsDouble(PyTuple_GET_ITEM(x_i_py, 1));
- }
- else {
- N = PyInt_AsLong(x_i_py);
- if (N==-1 && PyErr_Occurred()) goto fail;
- dx = 1.0;
- }
- equal = 1;
- }
- N -= 1;
-
- if (N < 2) {
- PyErr_Format(PyExc_ValueError, "too few samples (%d)", N);
- return NULL;
- }
- /* create output matrix */
- dims[0] = N-1;
- dims[1] = N+k;
- BB = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);
- if (BB == NULL) goto fail;
-
- t = (double *)malloc(sizeof(double)*(N+2*k-1));
- if (t==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
-
- /* Create work array to hold computed non-zero values for
- the spline for a value of x.
- */
- h = (double *)malloc(sizeof(double)*(2*k+1));
- if (h==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
-
- if (equal) { /* points equally spaced by 1 */
- /* we run deBoor's full derivative algorithm twice, subtract the results
- offset by one and then copy the result one time with artificially created knots
- Then, we keep copying the result to every row */
-
- /* Create knots at equally-spaced locations from -(K-1) to N+K-1 */
- double *tmp, factor;
- int numbytes;
- numbytes = (k+2)*sizeof(double);
- tmp = malloc(numbytes);
- if (tmp==NULL) {
- PyErr_NoMemory();
- goto fail;
- }
- ptr = t;
- for (i=-k+1; i<N+k; i++) *ptr++ = i;
- j = k-1;
- _deBoor_D(t, 0, k, k-1, k, h);
- ptr = tmp;
- for (m=0; m<=k; m++) *ptr++ = -h[m];
- _deBoor_D(t, 0, k, k, k, h);
- ptr = tmp+1;
- for (m=0; m<=k; m++) *ptr++ += h[m];
- if (dx != 1.0) {
- factor = pow(dx, (double)k);
- for (m=0; m<(k+2); m++) {
- tmp[m] /= factor;
- }
- }
- ptr = PyArray_DATA(BB);
- for (i=0; i<(N-1); i++) {
- memcpy(ptr, tmp, numbytes);
- ptr += (N+k+1);
- }
- free(tmp);
- goto finish;
- }
-
- /* Not-equally spaced */
- x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
- if (x_i == NULL) return NULL;
-
- /* create dummy knot array with new knots inserted at the end
- selected as mirror symmetric versions of the old knots
- */
- x0 = *((double *)PyArray_DATA(x_i));
- xN = *((double *)PyArray_DATA(x_i) + N);
- for (i=0; i<k-1; i++) { /* fill in ends if k > 1*/
- t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i,k-1-i)));
- t[k+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i,N-1-i)));
- }
- ptr = t + (k-1);
- for (i=0; i<=N; i++) {
- *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)));
- }
-
-
- /* Determine the K+1 non-zero values of the discontinuity jump matrix
- and place them in the correct location in the matrix for each row
- (along the diagonals).
-
- The matrix is
-
- J_{ij} = b^{(k)}_{j,k}(x^{+}_i) - b^{(k)}_{j,k}(x^{-}_i)
-
- */
- ptr = PyArray_DATA(BB);
- dptr = ptr;
- for (i=0,j=k-1; i<N-1; i++,j++) {
- _deBoor_D(t, 0, k, j, k, h);
- /* We need to copy over but negate the terms */
- for (m=0; m<=k; m++) *ptr++ = -h[m];
- /* If we are past the first row, then we need to also add the current
- values result to the previous row */
- if (i>0) {
- for (m=0; m<=k; m++) *dptr++ += h[m];
- }
- /* store location of last start position plus one.*/
- dptr = ptr - k;
- ptr += N; /* advance to next row shifted over one */
- }
- /* We need to finish the result for the last row. */
- _deBoor_D(t, 0, k, j, k, h);
- for (m=0; m<=k; m++) *dptr++ += h[m];
-
- finish:
- Py_XDECREF(x_i);
- free(t);
- free(h);
- return (PyObject *)BB;
-
- fail:
- Py_XDECREF(x_i);
- Py_XDECREF(BB);
- if (t != NULL) free(t);
- if (h != NULL) free(h);
- return NULL;
-}
-
-
-
Deleted: trunk/scipy/interpolate/_fitpackmodule.c
===================================================================
--- trunk/scipy/interpolate/_fitpackmodule.c 2008-09-04 16:47:55 UTC (rev 4686)
+++ trunk/scipy/interpolate/_fitpackmodule.c 2008-09-04 19:41:11 UTC (rev 4687)
@@ -1,36 +0,0 @@
-/*
- Multipack project.
- This file is generated by setmodules.py. Do not modify it.
- */
-#include "multipack.h"
-static PyObject *fitpack_error;
-#include "__fitpack.h"
-static struct PyMethodDef fitpack_module_methods[] = {
-{"_curfit", fitpack_curfit, METH_VARARGS, doc_curfit},
-{"_spl_", fitpack_spl_, METH_VARARGS, doc_spl_},
-{"_splint", fitpack_splint, METH_VARARGS, doc_splint},
-{"_sproot", fitpack_sproot, METH_VARARGS, doc_sproot},
-{"_spalde", fitpack_spalde, METH_VARARGS, doc_spalde},
-{"_parcur", fitpack_parcur, METH_VARARGS, doc_parcur},
-{"_surfit", fitpack_surfit, METH_VARARGS, doc_surfit},
-{"_bispev", fitpack_bispev, METH_VARARGS, doc_bispev},
-{"_insert", fitpack_insert, METH_VARARGS, doc_insert},
-{"_bspleval", _bspleval, METH_VARARGS, doc_bspleval},
-{"_bsplmat", _bsplmat, METH_VARARGS, doc_bsplmat},
-{"_bspldismat", _bspldismat, METH_VARARGS, doc_bspldismat},
-{NULL, NULL, 0, NULL}
-};
-PyMODINIT_FUNC init_fitpack(void) {
- PyObject *m, *d, *s;
- m = Py_InitModule("_fitpack", fitpack_module_methods);
- import_array();
- d = PyModule_GetDict(m);
-
- s = PyString_FromString(" 1.7 ");
- PyDict_SetItemString(d, "__version__", s);
- fitpack_error = PyErr_NewException ("fitpack.error", NULL, NULL);
- Py_DECREF(s);
- if (PyErr_Occurred())
- Py_FatalError("can't initialize module fitpack");
-}
-
Deleted: trunk/scipy/interpolate/fitpack.pyf
===================================================================
--- trunk/scipy/interpolate/fitpack.pyf 2008-09-04 16:47:55 UTC (rev 4686)
+++ trunk/scipy/interpolate/fitpack.pyf 2008-09-04 19:41:11 UTC (rev 4687)
@@ -1,479 +0,0 @@
-! -*- f90 -*-
-! Author: Pearu Peterson <pearu at cens.ioc.ee>
-!
-python module dfitpack ! in
-
- usercode '''
-
-static double dmax(double* seq,int len) {
- double val;
- int i;
- if (len<1)
- return -1e308;
- val = seq[0];
- for(i=1;i<len;++i)
- if (seq[i]>val) val = seq[i];
- return val;
-}
-static double dmin(double* seq,int len) {
- double val;
- int i;
- if (len<1)
- return 1e308;
- val = seq[0];
- for(i=1;i<len;++i)
- if (seq[i]<val) val = seq[i];
- return val;
-}
-static double calc_b(double* x,int m,double* tx,int nx) {
- double val1 = dmin(x,m);
- double val2 = dmin(tx,nx);
- if (val2>val1) return val1;
- val1 = dmax(tx,nx);
- return val2 - (val1-val2)/nx;
-}
-static double calc_e(double* x,int m,double* tx,int nx) {
- double val1 = dmax(x,m);
- double val2 = dmax(tx,nx);
- if (val2<val1) return val1;
- val1 = dmin(tx,nx);
- return val2 + (val2-val1)/nx;
-}
-static int imax(int i1,int i2) {
- return MAX(i1,i2);
-}
-
-static int calc_surfit_lwrk1(int m, int kx, int ky, int nxest, int nyest) {
- int u = nxest-kx-1;
- int v = nyest-ky-1;
- int km = MAX(kx,ky)+1;
- int ne = MAX(nxest,nyest);
- int bx = kx*v+ky+1;
- int by = ky*u+kx+1;
- int b1,b2;
- if (bx<=by) {b1=bx;b2=bx+v-ky;}
- else {b1=by;b2=by+u-kx;}
- return u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1;
-}
-static int calc_surfit_lwrk2(int m, int kx, int ky, int nxest, int nyest) {
- int u = nxest-kx-1;
- int v = nyest-ky-1;
- int bx = kx*v+ky+1;
- int by = ky*u+kx+1;
- int b2 = (bx<=by?bx+v-ky:by+u-kx);
- return u*v*(b2+1)+b2;
-}
-
-static int calc_regrid_lwrk(int mx, int my, int kx, int ky,
- int nxest, int nyest) {
- int u = MAX(my, nxest);
- return 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+my*(ky+1)+u;
-}
-'''
-
- interface
-
- !!!!!!!!!! Univariate spline !!!!!!!!!!!
-
- subroutine splev(t,n,c,k,x,y,m,ier)
- ! y = splev(t,c,k,x)
- real*8 dimension(n),intent(in) :: t
- integer intent(hide),depend(t) :: n=len(t)
- real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
- integer :: k
- real*8 dimension(m),intent(in) :: x
- real*8 dimension(m),depend(m),intent(out) :: y
- integer intent(hide),depend(x) :: m=len(x)
- integer intent(hide) :: ier
- end subroutine splev
-
- subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
- ! dy = splder(t,c,k,x,[nu])
- real*8 dimension(n) :: t
- integer depend(t),intent(hide) :: n=len(t)
- real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
- integer :: k
- integer depend(k),check(0<=nu && nu<=k) :: nu = 1
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),intent(out) :: y
- integer depend(x),intent(hide) :: m=len(x)
- real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
- integer intent(hide) :: ier
- end subroutine splder
-
- function splint(t,n,c,k,a,b,wrk)
- ! iy = splint(t,c,k,a,b)
- real*8 dimension(n),intent(in) :: t
- integer intent(hide),depend(t) :: n=len(t)
- real*8 dimension(n),depend(n),check(len(c)==n) :: c
- integer intent(in) :: k
- real*8 intent(in) :: a
- real*8 intent(in) :: b
- real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
- real*8 :: splint
- end function splint
-
- subroutine sproot(t,n,c,zero,mest,m,ier)
- ! zero,m,ier = sproot(t,c[,mest])
- real*8 dimension(n),intent(in) :: t
- integer intent(hide),depend(t),check(n>=8) :: n=len(t)
- real*8 dimension(n),depend(n),check(len(c)==n) :: c
- real*8 dimension(mest),intent(out),depend(mest) :: zero
- integer optional,intent(in),depend(n) :: mest=3*(n-7)
- integer intent(out) :: m
- integer intent(out) :: ier
- end subroutine sproot
-
- subroutine spalde(t,n,c,k,x,d,ier)
- ! d,ier = spalde(t,c,k,x)
-
- callprotoargument double*,int*,double*,int*,double*,double*,int*
- callstatement {int k1=k+1; (*f2py_func)(t,&n,c,&k1,&x,d,&ier); }
-
- real*8 dimension(n) :: t
- integer intent(hide),depend(t) :: n=len(t)
- real*8 dimension(n),depend(n),check(len(c)==n) :: c
- integer intent(in) :: k
- real*8 intent(in) :: x
- real*8 dimension(k+1),intent(out),depend(k) :: d
- integer intent(out) :: ier
- end subroutine spalde
-
- subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
- ! in curfit.f
- integer :: iopt
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(w)==m) :: w
- real*8 optional,depend(x),check(xb<=x[0]) :: xb = x[0]
- real*8 optional,depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
- integer optional,check(1<=k && k <=5),intent(in) :: k=3
- real*8 optional,check(s>=0.0) :: s = 0.0
- integer intent(hide),depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest),intent(inout) :: t
- real*8 dimension(n),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest),intent(inout) :: iwrk
- integer intent(out) :: ier
- end subroutine curfit
-
- subroutine percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
- ! in percur.f
- integer :: iopt
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(w)==m) :: w
- integer optional,check(1<=k && k <=5),intent(in) :: k=3
- real*8 optional,check(s>=0.0) :: s = 0.0
- integer intent(hide),depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest),intent(inout) :: t
- real*8 dimension(n),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest),intent(inout) :: iwrk
- integer intent(out) :: ier
- end subroutine percur
-
-
- subroutine parcur(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,c,fp,wrk,lwrk,iwrk,ier)
- ! in parcur.f
- integer check(iopt>=-1 && iopt <= 1):: iopt
- integer check(ipar == 1 || ipar == 0) :: ipar
- integer check(idim > 0 && idim < 11) :: idim
- integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
- real*8 dimension(m), intent(inout) :: u
- integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
- real*8 dimension(mx) :: x
- real*8 dimension(m) :: w
- real*8 :: ub
- real*8 :: ue
- integer optional, check(1<=k && k<=5) :: k=3.0
- real*8 optional, check(s>=0.0) :: s = 0.0
- integer intent(hide), depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest), intent(inout) :: t
- integer intent(hide), depend(c,nest,idim), check(nc>=idim*nest) :: nc=len(c)
- real*8 dimension(nc), intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk), intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest), intent(inout) :: iwrk
- integer intent(out) :: ier
- end subroutine parcur
-
-
- subroutine fpcurf0(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
- ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurf0(x,y,k,[w,xb,xe,s,nest])
-
- fortranname fpcurf
- callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
- callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
- integer intent(hide) :: iopt = 0
- real*8 dimension(m),intent(in,out) :: x
- real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
- real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
- real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
- integer check(1<=k && k<=5),intent(in,out) :: k
- real*8 check(s>=0.0),depend(m),intent(in,out) :: s = m
- integer intent(in),depend(m,s,k,k1),check(nest>=2*k1) :: nest = (s==0.0?m+k+1:MAX(m/2,2*k1))
- real*8 intent(hide) :: tol = 0.001
- integer intent(hide) :: maxit = 20
- integer intent(hide),depend(k) :: k1=k+1
- integer intent(hide),depend(k) :: k2=k+2
- integer intent(out) :: n
- real*8 dimension(nest),intent(out),depend(nest) :: t
- real*8 dimension(nest),depend(nest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(nest),depend(nest),intent(out,cache) :: fpint
- real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
- integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
- integer intent(out) :: ier
- end subroutine fpcurf0
-
- subroutine fpcurf1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
- ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurf1(x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier)
-
- fortranname fpcurf
- callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
- callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
- integer intent(hide) :: iopt = 1
- real*8 dimension(m),intent(in,out,overwrite) :: x
- real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out,overwrite) :: y
- real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out,overwrite) :: w
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 intent(in,out) :: xb
- real*8 intent(in,out) :: xe
- integer check(1<=k && k<=5),intent(in,out) :: k
- real*8 check(s>=0.0),intent(in,out) :: s
- integer intent(hide),depend(t) :: nest = len(t)
- real*8 intent(hide) :: tol = 0.001
- integer intent(hide) :: maxit = 20
- integer intent(hide),depend(k) :: k1=k+1
- integer intent(hide),depend(k) :: k2=k+2
- integer intent(in,out) :: n
- real*8 dimension(nest),intent(in,out,overwrite) :: t
- real*8 dimension(nest),depend(nest),check(len(c)==nest),intent(in,out,overwrite) :: c
- real*8 intent(in,out) :: fp
- real*8 dimension(nest),depend(nest),check(len(fpint)==nest),intent(in,out,cache,overwrite) :: fpint
- real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
- integer dimension(nest),depend(nest),check(len(nrdata)==nest),intent(in,out,cache,overwrite) :: nrdata
- integer intent(in,out) :: ier
- end subroutine fpcurf1
-
- subroutine fpcurfm1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
- ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurfm1(x,y,k,t,[w,xb,xe])
-
- fortranname fpcurf
- callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
- callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
- integer intent(hide) :: iopt = -1
- real*8 dimension(m),intent(in,out) :: x
- real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
- real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
- real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
- integer check(1<=k && k<=5),intent(in,out) :: k
- real*8 intent(out) :: s = -1
- integer intent(hide),depend(n) :: nest = n
- real*8 intent(hide) :: tol = 0.001
- integer intent(hide) :: maxit = 20
- integer intent(hide),depend(k) :: k1=k+1
- integer intent(hide),depend(k) :: k2=k+2
- integer intent(out),depend(t) :: n = len(t)
- real*8 dimension(n),intent(in,out,overwrite) :: t
- real*8 dimension(nest),depend(nest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(nest),depend(nest),intent(out,cache) :: fpint
- real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
- integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
- integer intent(out) :: ier
- end subroutine fpcurfm1
-
- !!!!!!!!!! Bivariate spline !!!!!!!!!!!
-
- subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,iwrk,kwrk,ier)
- ! z,ier = bispev(tx,ty,c,kx,ky,x,y)
- real*8 dimension(nx),intent(in) :: tx
- integer intent(hide),depend(tx) :: nx=len(tx)
- real*8 dimension(ny),intent(in) :: ty
- integer intent(hide),depend(ty) :: ny=len(ty)
- real*8 intent(in),dimension((nx-kx-1)*(ny-ky-1)),depend(nx,ny,kx,ky),&
- check(len(c)==(nx-kx-1)*(ny-ky-1)):: c
- integer :: kx
- integer :: ky
- real*8 intent(in),dimension(mx) :: x
- integer intent(hide),depend(x) :: mx=len(x)
- real*8 intent(in),dimension(my) :: y
- integer intent(hide),depend(y) :: my=len(y)
- real*8 dimension(mx,my),depend(mx,my),intent(out,c) :: z
- real*8 dimension(lwrk),depend(lwrk),intent(hide,cache) :: wrk
- integer intent(hide),depend(mx,kx,my,ky) :: lwrk=mx*(kx+1)+my*(ky+1)
- integer dimension(kwrk),depend(kwrk),intent(hide,cache) :: iwrk
- integer intent(hide),depend(mx,my) :: kwrk=mx+my
- integer intent(out) :: ier
- end subroutine bispev
-
- subroutine surfit_smth(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
- nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
- iwrk,kwrk,ier)
- ! nx,tx,ny,ty,c,fp,ier = surfit_smth(x,y,z,[w,xb,xe,yb,ye,kx,ky,s,eps,lwrk2])
-
- fortranname surfit
-
- integer intent(hide) :: iopt=0
- integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
- :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(z)==m) :: z
- real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
- real*8 optional,depend(x,m) :: xb=dmin(x,m)
- real*8 optional,depend(x,m) :: xe=dmax(x,m)
- real*8 optional,depend(y,m) :: yb=dmin(y,m)
- real*8 optional,depend(y,m) :: ye=dmax(y,m)
- integer check(1<=kx && kx<=5) :: kx = 3
- integer check(1<=ky && ky<=5) :: ky = 3
- real*8 optional,check(0.0<=s) :: s = m
- integer optional,depend(kx,m),check(nxest>=2*(kx+1)) &
- :: nxest = imax(kx+1+sqrt(m/2),2*(kx+1))
- integer optional,depend(ky,m),check(nyest>=2*(ky+1)) &
- :: nyest = imax(ky+1+sqrt(m/2),2*(ky+1))
- integer intent(hide),depend(nxest,nyest) :: nmax=MAX(nxest,nyest)
- real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
- integer intent(out) :: nx
- real*8 dimension(nmax),intent(out),depend(nmax) :: tx
- integer intent(out) :: ny
- real*8 dimension(nmax),intent(out),depend(nmax) :: ty
- real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
- depend(kx,ky,nxest,nyest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk1),intent(cache,out),depend(lwrk1) :: wrk1
- integer intent(hide),depend(m,kx,ky,nxest,nyest) &
- :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
- real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
- integer optional,intent(in),depend(kx,ky,nxest,nyest) &
- :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
- integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
- integer intent(hide),depend(m,nxest,nyest,kx,ky) &
- :: kwrk=m+(nxest-2*kx-1)*(nyest-2*ky-1)
- integer intent(out) :: ier
- end subroutine surfit_smth
-
- subroutine surfit_lsq(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
- nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
- iwrk,kwrk,ier)
- ! tx,ty,c,fp,ier = surfit_lsq(x,y,z,tx,ty,[w,xb,xe,yb,ye,kx,ky,eps,lwrk2])
-
- fortranname surfit
-
- integer intent(hide) :: iopt=-1
- integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
- :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(z)==m) :: z
- real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
- real*8 optional,depend(x,tx,m,nx) :: xb=calc_b(x,m,tx,nx)
- real*8 optional,depend(x,tx,m,nx) :: xe=calc_e(x,m,tx,nx)
- real*8 optional,depend(y,ty,m,ny) :: yb=calc_b(y,m,ty,ny)
- real*8 optional,depend(y,ty,m,ny) :: ye=calc_e(y,m,ty,ny)
- integer check(1<=kx && kx<=5) :: kx = 3
- integer check(1<=ky && ky<=5) :: ky = 3
- real*8 intent(hide) :: s = 0.0
- integer intent(hide),depend(nx) :: nxest = nx
- integer intent(hide),depend(ny) :: nyest = ny
- integer intent(hide),depend(nx,ny) :: nmax=MAX(nx,ny)
- real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
- integer intent(hide),depend(tx,kx),check(2*kx+2<=nx) :: nx = len(tx)
- real*8 dimension(nx),intent(in,out,overwrite) :: tx
- integer intent(hide),depend(ty,ky),check(2*ky+2<=ny) :: ny = len(ty)
- real*8 dimension(ny),intent(in,out,overwrite) :: ty
- real*8 dimension((nx-kx-1)*(ny-ky-1)),depend(kx,ky,nx,ny),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk1),intent(cache,hide),depend(lwrk1) :: wrk1
- integer intent(hide),depend(m,kx,ky,nxest,nyest) &
- :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
- real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
- integer optional,intent(in),depend(kx,ky,nxest,nyest) &
- :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
- integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
- integer intent(hide),depend(m,nx,ny,kx,ky) &
- :: kwrk=m+(nx-2*kx-1)*(ny-2*ky-1)
- integer intent(out) :: ier
- end subroutine surfit_lsq
-
- subroutine regrid_smth(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,&
- nxest,nyest,nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
- ! nx,tx,ny,ty,c,fp,ier = regrid_smth(x,y,z,[xb,xe,yb,ye,kx,ky,s])
-
- fortranname regrid
-
- integer intent(hide) :: iopt=0
- integer intent(hide),depend(x,kx),check(mx>kx) :: mx=len(x)
- real*8 dimension(mx) :: x
- integer intent(hide),depend(y,ky),check(my>ky) :: my=len(y)
- real*8 dimension(my) :: y
- real*8 dimension(mx*my),depend(mx,my),check(len(z)==mx*my) :: z
- real*8 optional,depend(x,mx) :: xb=dmin(x,mx)
- real*8 optional,depend(x,mx) :: xe=dmax(x,mx)
- real*8 optional,depend(y,my) :: yb=dmin(y,my)
- real*8 optional,depend(y,my) :: ye=dmax(y,my)
- integer optional,check(1<=kx && kx<=5) :: kx = 3
- integer optional,check(1<=ky && ky<=5) :: ky = 3
- real*8 optional,check(0.0<=s) :: s = 0.0
- integer intent(hide),depend(kx,mx),check(nxest>=2*(kx+1)) &
- :: nxest = mx+kx+1
- integer intent(hide),depend(ky,my),check(nyest>=2*(ky+1)) &
- :: nyest = my+ky+1
- integer intent(out) :: nx
- real*8 dimension(nxest),intent(out),depend(nxest) :: tx
- integer intent(out) :: ny
- real*8 dimension(nyest),intent(out),depend(nyest) :: ty
- real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
- depend(kx,ky,nxest,nyest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(cache,hide),depend(lwrk) :: wrk
- integer intent(hide),depend(mx,my,kx,ky,nxest,nyest) &
- :: lwrk=calc_regrid_lwrk(mx,my,kx,ky,nxest,nyest)
- integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
- integer intent(hide),depend(mx,my,nxest,nyest) &
- :: kwrk=3+mx+my+nxest+nyest
- integer intent(out) :: ier
- end subroutine regrid_smth
-
- function dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
- ! iy = dblint(tx,ty,c,kx,ky,xb,xe,yb,ye)
- real*8 dimension(nx),intent(in) :: tx
- integer intent(hide),depend(tx) :: nx=len(tx)
- real*8 dimension(ny),intent(in) :: ty
- integer intent(hide),depend(ty) :: ny=len(ty)
- real*8 intent(in),dimension((nx-kx-1)*(ny-ky-1)),depend(nx,ny,kx,ky),&
- check(len(c)==(nx-kx-1)*(ny-ky-1)):: c
- integer :: kx
- integer :: ky
- real*8 intent(in) :: xb
- real*8 intent(in) :: xe
- real*8 intent(in) :: yb
- real*8 intent(in) :: ye
- real*8 dimension(nx+ny-kx-ky-2),depend(nx,ny,kx,ky),intent(cache,hide) :: wrk
- real*8 :: dblint
- end function dblint
- end interface
-end python module dfitpack
-
Deleted: trunk/scipy/interpolate/multipack.h
===================================================================
--- trunk/scipy/interpolate/multipack.h 2008-09-04 16:47:55 UTC (rev 4686)
+++ trunk/scipy/interpolate/multipack.h 2008-09-04 19:41:11 UTC (rev 4687)
@@ -1,211 +0,0 @@
-/* MULTIPACK module by Travis Oliphant
-
-Copyright (c) 2002 Travis Oliphant all rights reserved
-Oliphant.Travis at altavista.net
-Permission to use, modify, and distribute this software is given under the
-terms of the SciPy (BSD style) license. See LICENSE.txt that came with
-this distribution for specifics.
-
-NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
-*/
-
-
-/* This extension module is a collection of wrapper functions around
-common FORTRAN code in the packages MINPACK, ODEPACK, and QUADPACK plus
-some differential algebraic equation solvers.
-
-The wrappers are meant to be nearly direct translations between the
-FORTAN code and Python. Some parameters like sizes do not need to be
-passed since they are available from the objects.
-
-It is anticipated that a pure Python module be written to call these lower
-level routines and make a simpler user interface. All of the routines define
-default values for little-used parameters so that even the raw routines are
-quite useful without a separate wrapper.
-
-FORTRAN Outputs that are not either an error indicator or the sought-after
-results are placed in a dictionary and returned as an optional member of
-the result tuple when the full_output argument is non-zero.
-*/
-
-#include "Python.h"
-#include "numpy/arrayobject.h"
-
-#define PYERR(errobj,message) {PyErr_SetString(errobj,message); goto fail;}
-#define PYERR2(errobj,message) {PyErr_Print(); PyErr_SetString(errobj, message); goto fail;}
-#define ISCONTIGUOUS(m) ((m)->flags & CONTIGUOUS)
-
-#define STORE_VARS() PyObject *store_multipack_globals[4]; int store_multipack_globals3;
-
-#define INIT_FUNC(fun,arg,errobj) { /* Get extra arguments or set to zero length tuple */ \
- store_multipack_globals[0] = multipack_python_function; \
- store_multipack_globals[1] = multipack_extra_arguments; \
- if (arg == NULL) { \
- if ((arg = PyTuple_New(0)) == NULL) goto fail; \
- } \
- else \
- Py_INCREF(arg); /* We decrement on exit. */ \
- if (!PyTuple_Check(arg)) \
- PYERR(errobj,"Extra Arguments must be in a tuple"); \
- /* Set up callback functions */ \
- if (!PyCallable_Check(fun)) \
- PYERR(errobj,"First argument must be a callable function."); \
- multipack_python_function = fun; \
- multipack_extra_arguments = arg; }
-
-#define INIT_JAC_FUNC(fun,Dfun,arg,col_deriv,errobj) { \
- store_multipack_globals[0] = multipack_python_function; \
- store_multipack_globals[1] = multipack_extra_arguments; \
- store_multipack_globals[2] = multipack_python_jacobian; \
- store_multipack_globals3 = multipack_jac_transpose; \
- if (arg == NULL) { \
- if ((arg = PyTuple_New(0)) == NULL) goto fail; \
- } \
- else \
- Py_INCREF(arg); /* We decrement on exit. */ \
- if (!PyTuple_Check(arg)) \
- PYERR(errobj,"Extra Arguments must be in a tuple"); \
- /* Set up callback functions */ \
- if (!PyCallable_Check(fun) || (Dfun != Py_None && !PyCallable_Check(Dfun))) \
- PYERR(errobj,"The function and its Jacobian must be callable functions."); \
- multipack_python_function = fun; \
- multipack_extra_arguments = arg; \
- multipack_python_jacobian = Dfun; \
- multipack_jac_transpose = !(col_deriv);}
-
-#define RESTORE_JAC_FUNC() multipack_python_function = store_multipack_globals[0]; \
- multipack_extra_arguments = store_multipack_globals[1]; \
- multipack_python_jacobian = store_multipack_globals[2]; \
- multipack_jac_transpose = store_multipack_globals3;
-
-#define RESTORE_FUNC() multipack_python_function = store_multipack_globals[0]; \
- multipack_extra_arguments = store_multipack_globals[1];
-
-#define SET_DIAG(ap_diag,o_diag,mode) { /* Set the diag vector from input */ \
- if (o_diag == NULL || o_diag == Py_None) { \
- ap_diag = (PyArrayObject *)PyArray_FromDims(1,&n,PyArray_DOUBLE); \
- if (ap_diag == NULL) goto fail; \
- diag = (double *)ap_diag -> data; \
- mode = 1; \
- } \
- else { \
- ap_diag = (PyArrayObject *)PyArray_ContiguousFromObject(o_diag, PyArray_DOUBLE, 1, 1); \
- if (ap_diag == NULL) goto fail; \
- diag = (double *)ap_diag -> data; \
- mode = 2; } }
-
-#define MATRIXC2F(jac,data,n,m) {double *p1=(double *)(jac), *p2, *p3=(double *)(data);\
-int i,j;\
-for (j=0;j<(m);p3++,j++) \
- for (p2=p3,i=0;i<(n);p2+=(m),i++,p1++) \
- *p1 = *p2; }
-/*
-static PyObject *multipack_python_function=NULL;
-static PyObject *multipack_python_jacobian=NULL;
-static PyObject *multipack_extra_arguments=NULL;
-static int multipack_jac_transpose=1;
-*/
-
-static PyArrayObject * my_make_numpy_array(PyObject *y0, int type, int mindim, int maxdim)
- /* This is just like PyArray_ContiguousFromObject except it handles
- * single numeric datatypes as 1-element, rank-1 arrays instead of as
- * scalars.
- */
-{
- PyArrayObject *new_array;
- PyObject *tmpobj;
-
- Py_INCREF(y0);
-
- if (PyInt_Check(y0) || PyFloat_Check(y0)) {
- tmpobj = PyList_New(1);
- PyList_SET_ITEM(tmpobj, 0, y0); /* reference now belongs to tmpobj */
- }
- else
- tmpobj = y0;
-
- new_array = (PyArrayObject *)PyArray_ContiguousFromObject(tmpobj, type, mindim, maxdim);
-
- Py_DECREF(tmpobj);
- return new_array;
-}
-
-static PyObject *call_python_function(PyObject *func, int n, double *x, PyObject *args, int dim, PyObject *error_obj)
-{
- /*
- This is a generic function to call a python function that takes a 1-D
- sequence as a first argument and optional extra_arguments (should be a
- zero-length tuple if none desired). The result of the function is
- returned in a multiarray object.
- -- build sequence object from values in x.
- -- add extra arguments (if any) to an argument list.
- -- call Python callable object
- -- check if error occurred:
- if so return NULL
- -- if no error, place result of Python code into multiarray object.
- */
-
- PyArrayObject *sequence = NULL;
- PyObject *arglist = NULL, *tmpobj = NULL;
- PyObject *arg1 = NULL, *str1 = NULL;
- PyObject *result = NULL;
- PyArrayObject *result_array = NULL;
-
- /* Build sequence argument from inputs */
- sequence = (PyArrayObject *)PyArray_FromDimsAndData(1, &n, PyArray_DOUBLE, (char *)x);
- if (sequence == NULL) PYERR2(error_obj,"Internal failure to make an array of doubles out of first\n argument to function call.");
-
- /* Build argument list */
- if ((arg1 = PyTuple_New(1)) == NULL) {
- Py_DECREF(sequence);
- return NULL;
- }
- PyTuple_SET_ITEM(arg1, 0, (PyObject *)sequence);
- /* arg1 now owns sequence reference */
- if ((arglist = PySequence_Concat( arg1, args)) == NULL)
- PYERR2(error_obj,"Internal error constructing argument list.");
-
- Py_DECREF(arg1); /* arglist has a reference to sequence, now. */
-
-
- /* Call function object --- variable passed to routine. Extra
- arguments are in another passed variable.
- */
- if ((result = PyEval_CallObject(func, arglist))==NULL) {
- PyErr_Print();
- tmpobj = PyObject_GetAttrString(func, "func_name");
- if (tmpobj == NULL) goto fail;
- str1 = PyString_FromString("Error occured while calling the Python function named ");
- if (str1 == NULL) { Py_DECREF(tmpobj); goto fail;}
- PyString_ConcatAndDel(&str1, tmpobj);
- PyErr_SetString(error_obj, PyString_AsString(str1));
- Py_DECREF(str1);
- goto fail;
- }
-
- if ((result_array = (PyArrayObject *)PyArray_ContiguousFromObject(result, PyArray_DOUBLE, dim-1, dim))==NULL)
- PYERR2(error_obj,"Result from function call is not a proper array of floats.");
-
- Py_DECREF(result);
- Py_DECREF(arglist);
- return (PyObject *)result_array;
-
- fail:
- Py_XDECREF(arglist);
- Py_XDECREF(result);
- Py_XDECREF(arg1);
- return NULL;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
Modified: trunk/scipy/interpolate/setup.py
===================================================================
--- trunk/scipy/interpolate/setup.py 2008-09-04 16:47:55 UTC (rev 4686)
+++ trunk/scipy/interpolate/setup.py 2008-09-04 19:41:11 UTC (rev 4687)
@@ -12,12 +12,12 @@
)
config.add_extension('_fitpack',
- sources=['_fitpackmodule.c'],
+ sources=['src/_fitpackmodule.c'],
libraries=['fitpack'],
)
config.add_extension('dfitpack',
- sources=['fitpack.pyf'],
+ sources=['src/fitpack.pyf'],
libraries=['fitpack'],
)
Copied: trunk/scipy/interpolate/src/__fitpack.h (from rev 4684, trunk/scipy/interpolate/__fitpack.h)
Copied: trunk/scipy/interpolate/src/_fitpackmodule.c (from rev 4684, trunk/scipy/interpolate/_fitpackmodule.c)
Copied: trunk/scipy/interpolate/src/fitpack.pyf (from rev 4684, trunk/scipy/interpolate/fitpack.pyf)
Copied: trunk/scipy/interpolate/src/multipack.h (from rev 4684, trunk/scipy/interpolate/multipack.h)
More information about the Scipy-svn
mailing list