Bug Summary

File:/tmp/pyrefcon/scipy/scipy/interpolate/src/_fitpackmodule.c
Warning:line 1554, column 9
PyObject ownership leak with reference count of 1

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -disable-llvm-verifier -discard-value-names -main-file-name _fitpackmodule.c -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -analyzer-output=html -analyzer-checker=python -analyzer-disable-checker=deadcode -analyzer-config prune-paths=true,suppress-c++-stdlib=true,suppress-null-return-paths=false,crosscheck-with-z3=true,model-path=/opt/pyrefcon/lib/pyrefcon/models/models -analyzer-config experimental-enable-naive-ctu-analysis=true,ctu-dir=/tmp/pyrefcon/scipy/csa-scan,ctu-index-name=/tmp/pyrefcon/scipy/csa-scan/externalDefMap.txt,ctu-invocation-list=/tmp/pyrefcon/scipy/csa-scan/invocations.yaml,display-ctu-progress=false -setup-static-analyzer -analyzer-config-compatibility-mode=true -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -tune-cpu generic -debug-info-kind=limited -dwarf-version=4 -debugger-tuning=gdb -fcoverage-compilation-dir=/tmp/pyrefcon/scipy -resource-dir /opt/pyrefcon/lib/clang/13.0.0 -isystem /opt/pyrefcon/lib/pyrefcon/models/python3.8 -D NDEBUG -D _FORTIFY_SOURCE=2 -D NPY_NO_DEPRECATED_API=NPY_1_9_API_VERSION -I /usr/lib/python3/dist-packages/numpy/core/include -internal-isystem /opt/pyrefcon/lib/clang/13.0.0/include -internal-isystem /usr/local/include -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/10/../../../../x86_64-linux-gnu/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -Wno-unused-result -Wsign-compare -Wall -Wformat -Werror=format-security -Wformat -Werror=format-security -Wdate-time -fdebug-compilation-dir=/tmp/pyrefcon/scipy -ferror-limit 19 -fwrapv -pthread -stack-protector 2 -fgnuc-version=4.2.1 -vectorize-loops -vectorize-slp -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /tmp/pyrefcon/scipy/csa-scan/reports -x c scipy/interpolate/src/_fitpackmodule.c

scipy/interpolate/src/_fitpackmodule.c

→
1/*
2 Multipack project.
3 This file is generated by setmodules.py. Do not modify it.
4 */
5
6#include <Python.h>
7#include "numpy/arrayobject.h"
8
9#define PyInt_AsLongPyLong_AsLong PyLong_AsLong
10
11static PyObject *fitpack_error;
12#include "__fitpack.h"
13
14#ifdef HAVE_ILP64
15
16#define F_INTint npy_int64
17#define F_INT_NPYNPY_INT NPY_INT64NPY_LONG
18
19#if NPY_BITSOF_SHORT(2 * 8) == 64
20#define F_INT_PYFMT"i" "h"
21#elif NPY_BITSOF_INT(4 * 8) == 64
22#define F_INT_PYFMT"i" "i"
23#elif NPY_BITSOF_LONG(8 * 8) == 64
24#define F_INT_PYFMT"i" "l"
25#elif NPY_BITSOF_LONGLONG(8 * 8) == 64
26#define F_INT_PYFMT"i" "L"
27#else
28#error No compatible 64-bit integer size. \
29 Please contact NumPy maintainers and give detailed information about your \
30 compiler and platform, or set NPY_USE_BLAS64_=0
31#endif
32
33#else
34
35#define F_INTint int
36#define F_INT_NPYNPY_INT NPY_INT
37#define F_INT_PYFMT"i" "i"
38
39#endif
40
41
42/*
43 * Functions moved verbatim from __fitpack.h
44 */
45
46
47/*
48 * Python-C wrapper of FITPACK (by P. Dierckx) (in netlib known as dierckx)
49 * Author: Pearu Peterson <pearu@ioc.ee>
50 * June 1.-4., 1999
51 * June 7. 1999
52 * $Revision$
53 * $Date$
54 */
55
56/* module_methods:
57 * {"_curfit", fitpack_curfit, METH_VARARGS, doc_curfit},
58 * {"_spl_", fitpack_spl_, METH_VARARGS, doc_spl_},
59 * {"_splint", fitpack_splint, METH_VARARGS, doc_splint},
60 * {"_sproot", fitpack_sproot, METH_VARARGS, doc_sproot},
61 * {"_spalde", fitpack_spalde, METH_VARARGS, doc_spalde},
62 * {"_parcur", fitpack_parcur, METH_VARARGS, doc_parcur},
63 * {"_surfit", fitpack_surfit, METH_VARARGS, doc_surfit},
64 * {"_bispev", fitpack_bispev, METH_VARARGS, doc_bispev},
65 * {"_insert", fitpack_insert, METH_VARARGS, doc_insert},
66 */
67
68/* link libraries: (one item per line)
69 ddierckx
70 */
71/* python files: (to be imported to Multipack.py)
72 fitpack.py
73 */
74
75#if defined(UPPERCASE_FORTRAN)
76 #if defined(NO_APPEND_FORTRAN)
77 /* nothing to do */
78 #else
79 #define CURFITcurfit_ CURFIT_
80 #define PERCURpercur_ PERCUR_
81 #define SPALDEspalde_ SPALDE_
82 #define SPLDERsplder_ SPLDER_
83 #define SPLEVsplev_ SPLEV_
84 #define SPLINTsplint_ SPLINT_
85 #define SPROOTsproot_ SPROOT_
86 #define PARCURparcur_ PARCUR_
87 #define CLOCURclocur_ CLOCUR_
88 #define SURFITsurfit_ SURFIT_
89 #define BISPEVbispev_ BISPEV_
90 #define PARDERparder_ PARDER_
91 #define INSERTinsert_ INSERT_
92 #endif
93#else
94 #if defined(NO_APPEND_FORTRAN)
95 #define CURFITcurfit_ curfit
96 #define PERCURpercur_ percur
97 #define SPALDEspalde_ spalde
98 #define SPLDERsplder_ splder
99 #define SPLEVsplev_ splev
100 #define SPLINTsplint_ splint
101 #define SPROOTsproot_ sproot
102 #define PARCURparcur_ parcur
103 #define CLOCURclocur_ clocur
104 #define SURFITsurfit_ surfit
105 #define BISPEVbispev_ bispev
106 #define PARDERparder_ parder
107 #define INSERTinsert_ insert
108 #else
109 #define CURFITcurfit_ curfit_
110 #define PERCURpercur_ percur_
111 #define SPALDEspalde_ spalde_
112 #define SPLDERsplder_ splder_
113 #define SPLEVsplev_ splev_
114 #define SPLINTsplint_ splint_
115 #define SPROOTsproot_ sproot_
116 #define PARCURparcur_ parcur_
117 #define CLOCURclocur_ clocur_
118 #define SURFITsurfit_ surfit_
119 #define BISPEVbispev_ bispev_
120 #define PARDERparder_ parder_
121 #define INSERTinsert_ insert_
122 #endif
123#endif
124
125void CURFITcurfit_(F_INTint*,F_INTint*,double*,double*,double*,double*,
126 double*,F_INTint*,double*,F_INTint*,F_INTint*,double*,double*,
127 double*,double*,F_INTint*,F_INTint*,F_INTint*);
128void PERCURpercur_(F_INTint*,F_INTint*,double*,double*,double*,F_INTint*,
129 double*,F_INTint*,F_INTint*,double*,double*,double*,
130 double*,F_INTint*,F_INTint*,F_INTint*);
131void SPALDEspalde_(double*,F_INTint*,double*,F_INTint*,double*,double*,F_INTint*);
132void SPLDERsplder_(double*,F_INTint*,double*,F_INTint*,F_INTint*,double*,
133 double*,F_INTint*,F_INTint*,double*,F_INTint*);
134void SPLEVsplev_(double*,F_INTint*,double*,F_INTint*,double*,double*,F_INTint*,F_INTint*,F_INTint*);
135double SPLINTsplint_(double*,F_INTint*,double*,F_INTint*,double*,double*,double*);
136void SPROOTsproot_(double*,F_INTint*,double*,double*,F_INTint*,F_INTint*,F_INTint*);
137void PARCURparcur_(F_INTint*,F_INTint*,F_INTint*,F_INTint*,double*,F_INTint*,double*,
138 double*,double*,double*,F_INTint*,double*,F_INTint*,F_INTint*,
139 double*,F_INTint*,double*,double*,double*,F_INTint*,F_INTint*,F_INTint*);
140void CLOCURclocur_(F_INTint*,F_INTint*,F_INTint*,F_INTint*,double*,F_INTint*,double*,
141 double*,F_INTint*,double*,F_INTint*,F_INTint*,double*,F_INTint*,
142 double*,double*,double*,F_INTint*,F_INTint*,F_INTint*);
143void SURFITsurfit_(F_INTint*,F_INTint*,double*,double*,double*,double*,
144 double*,double*,double*,double*,F_INTint*,F_INTint*,double*,
145 F_INTint*,F_INTint*,F_INTint*,double*,F_INTint*,double*,F_INTint*,double*,
146 double*,double*,double*,F_INTint*,double*,F_INTint*,F_INTint*,F_INTint*,F_INTint*);
147void BISPEVbispev_(double*,F_INTint*,double*,F_INTint*,double*,F_INTint*,F_INTint*,
148 double*,F_INTint*,double*,F_INTint*,double*,double*,F_INTint*,
149 F_INTint*,F_INTint*,F_INTint*);
150void PARDERparder_(double*,F_INTint*,double*,F_INTint*,double*,F_INTint*,F_INTint*,
151 F_INTint*,F_INTint*,double*,F_INTint*,double*,F_INTint*,double*,
152 double*,F_INTint*,F_INTint*,F_INTint*,F_INTint*);
153void INSERTinsert_(F_INTint*,double*,F_INTint*,double*,F_INTint*,double*,double*,
154 F_INTint*,double*,F_INTint*,F_INTint*);
155
156/* Note that curev, cualde need no interface. */
157
158static char doc_bispev[] = " [z,ier] = _bispev(tx,ty,c,kx,ky,x,y,nux,nuy)";
159static PyObject *
160fitpack_bispev(PyObject *dummy, PyObject *args)
161{
162 F_INTint nx, ny, kx, ky, mx, my, lwrk, *iwrk, kwrk, ier, lwa, nux, nuy;
163 npy_intp mxy;
164 double *tx, *ty, *c, *x, *y, *z, *wrk, *wa = NULL((void*)0);
165 PyArrayObject *ap_x = NULL((void*)0), *ap_y = NULL((void*)0), *ap_z = NULL((void*)0), *ap_tx = NULL((void*)0);
166 PyArrayObject *ap_ty = NULL((void*)0), *ap_c = NULL((void*)0);
167 PyObject *x_py = NULL((void*)0), *y_py = NULL((void*)0), *c_py = NULL((void*)0), *tx_py = NULL((void*)0), *ty_py = NULL((void*)0);
168
169 if (!PyArg_ParseTuple(args, ("OOO" F_INT_PYFMT"i" F_INT_PYFMT"i" "OO" F_INT_PYFMT"i" F_INT_PYFMT"i"),
170 &tx_py,&ty_py,&c_py,&kx,&ky,&x_py,&y_py,&nux,&nuy)) {
171 return NULL((void*)0);
172 }
173 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
174 ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(y_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
175 ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(c_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
176 ap_tx = (PyArrayObject *)PyArray_ContiguousFromObject(tx_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(tx_py, (*(PyArray_Descr * (*)(
int)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 |
0x0400))) | 0x0040, ((void*)0));
177 ap_ty = (PyArrayObject *)PyArray_ContiguousFromObject(ty_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(ty_py, (*(PyArray_Descr * (*)(
int)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 |
0x0400))) | 0x0040, ((void*)0));
178 if (ap_x == NULL((void*)0)
179 || ap_y == NULL((void*)0)
180 || ap_c == NULL((void*)0)
181 || ap_tx == NULL((void*)0)
182 || ap_ty == NULL((void*)0)) {
183 goto fail;
184 }
185 x = (double *) PyArray_DATA(ap_x);
186 y = (double *) PyArray_DATA(ap_y);
187 c = (double *) PyArray_DATA(ap_c);
188 tx = (double *) PyArray_DATA(ap_tx);
189 ty = (double *) PyArray_DATA(ap_ty);
190 nx = PyArray_DIMS(ap_tx)[0];
191 ny = PyArray_DIMS(ap_ty)[0];
192 mx = PyArray_DIMS(ap_x)[0];
193 my = PyArray_DIMS(ap_y)[0];
194 mxy = mx*my;
195 if (my != 0 && mxy/my != mx) {
196 /* Integer overflow */
197 PyErr_Format(PyExc_RuntimeError,
198 "Cannot produce output of size %dx%d (size too large)",
199 mx, my);
200 goto fail;
201 }
202 ap_z = (PyArrayObject *)PyArray_SimpleNew(1,&mxy,NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, &mxy, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
203 if (ap_z == NULL((void*)0)) {
204 goto fail;
205 }
206 z = (double *) PyArray_DATA(ap_z);
207 if (nux || nuy) {
208 lwrk = mx*(kx + 1 - nux) + my*(ky + 1 - nuy) + (nx - kx - 1)*(ny - ky - 1);
209 }
210 else {
211 lwrk = mx*(kx + 1) + my*(ky + 1);
212 }
213 kwrk = mx + my;
214 lwa = lwrk + kwrk;
215 if ((wa = malloc(lwa*sizeof(double))) == NULL((void*)0)) {
216 PyErr_NoMemory();
217 goto fail;
218 }
219 wrk = wa;
220 iwrk = (int *)(wrk + lwrk);
221 if (nux || nuy) {
222 PARDERparder_(tx, &nx, ty, &ny, c, &kx, &ky, &nux, &nuy, x, &mx, y, &my, z,
223 wrk, &lwrk, iwrk, &kwrk, &ier);
224 }
225 else {
226 BISPEVbispev_(tx, &nx, ty, &ny, c, &kx, &ky, x, &mx, y, &my, z, wrk, &lwrk,
227 iwrk, &kwrk, &ier);
228 }
229
230 free(wa);
231 Py_DECREF(ap_x)_Py_DECREF(((PyObject*)(ap_x)));
232 Py_DECREF(ap_y)_Py_DECREF(((PyObject*)(ap_y)));
233 Py_DECREF(ap_c)_Py_DECREF(((PyObject*)(ap_c)));
234 Py_DECREF(ap_tx)_Py_DECREF(((PyObject*)(ap_tx)));
235 Py_DECREF(ap_ty)_Py_DECREF(((PyObject*)(ap_ty)));
236 return Py_BuildValue("Ni",PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_z),ier);
237
238fail:
239 free(wa);
240 Py_XDECREF(ap_x)_Py_XDECREF(((PyObject*)(ap_x)));
241 Py_XDECREF(ap_y)_Py_XDECREF(((PyObject*)(ap_y)));
242 Py_XDECREF(ap_z)_Py_XDECREF(((PyObject*)(ap_z)));
243 Py_XDECREF(ap_c)_Py_XDECREF(((PyObject*)(ap_c)));
244 Py_XDECREF(ap_tx)_Py_XDECREF(((PyObject*)(ap_tx)));
245 Py_XDECREF(ap_ty)_Py_XDECREF(((PyObject*)(ap_ty)));
246 return NULL((void*)0);
247}
248
249static char doc_surfit[] = " [tx,ty,c,o] = _surfit(x, y, z, w, xb, xe, yb, ye,"\
250 " kx,ky,iopt,s,eps,tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)";
251static PyObject *
252fitpack_surfit(PyObject *dummy, PyObject *args)
253{
254 F_INTint iopt, m, kx, ky, nxest, nyest, lwrk1, lwrk2, *iwrk, kwrk, ier;
255 F_INTint lwa, nxo, nyo, i, lcest, nmax, nx, ny, lc;
256 npy_intp dims[1];
257 double *x, *y, *z, *w, xb, xe, yb, ye, s, *tx, *ty, *c, fp;
258 double *wrk1, *wrk2, *wa = NULL((void*)0), eps;
259 PyArrayObject *ap_x = NULL((void*)0), *ap_y = NULL((void*)0), *ap_z, *ap_w = NULL((void*)0);
260 PyArrayObject *ap_tx = NULL((void*)0),*ap_ty = NULL((void*)0),*ap_c = NULL((void*)0), *ap_wrk = NULL((void*)0);
261 PyObject *x_py = NULL((void*)0), *y_py = NULL((void*)0), *z_py = NULL((void*)0), *w_py = NULL((void*)0);
262 PyObject *tx_py = NULL((void*)0), *ty_py = NULL((void*)0), *wrk_py = NULL((void*)0);
263
264 nx = ny = ier = nxo = nyo = 0;
265 if (!PyArg_ParseTuple(args, ("OOOOdddd" F_INT_PYFMT"i" F_INT_PYFMT"i" F_INT_PYFMT"i"
266 "ddOO" F_INT_PYFMT"i" F_INT_PYFMT"i" "O"
267 F_INT_PYFMT"i" F_INT_PYFMT"i"),
268 &x_py, &y_py, &z_py, &w_py, &xb, &xe, &yb, &ye,
269 &kx, &ky, &iopt, &s, &eps, &tx_py, &ty_py, &nxest,
270 &nyest, &wrk_py, &lwrk1, &lwrk2)) {
271 return NULL((void*)0);
272 }
273 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
274 ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(y_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
275 ap_z = (PyArrayObject *)PyArray_ContiguousFromObject(z_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(z_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
276 ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(w_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
277 ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(wrk_py, (*(PyArray_Descr * (*)
(int)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100
| 0x0400))) | 0x0040, ((void*)0));
278 /*ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, F_INT_NPY, 0, 1);*/
279 if (ap_x == NULL((void*)0)
280 || ap_y == NULL((void*)0)
281 || ap_z == NULL((void*)0)
282 || ap_w == NULL((void*)0)
283 || ap_wrk == NULL((void*)0)) {
284 goto fail;
285 }
286 x = (double *) PyArray_DATA(ap_x);
287 y = (double *) PyArray_DATA(ap_y);
288 z = (double *) PyArray_DATA(ap_z);
289 w = (double *) PyArray_DATA(ap_w);
290 m = PyArray_DIMS(ap_x)[0];
291 nmax = nxest;
292 if (nmax < nyest) {
293 nmax = nyest;
294 }
295 lcest = (nxest - kx - 1)*(nyest - ky - 1);
296 kwrk = m + (nxest - 2*kx - 1)*(nyest - 2*ky - 1);
297 lwa = 2*nmax + lcest + lwrk1 + lwrk2 + kwrk;
298 if ((wa = malloc(lwa*sizeof(double))) == NULL((void*)0)) {
299 PyErr_NoMemory();
300 goto fail;
301 }
302 /*
303 * NOTE: the work arrays MUST be aligned on double boundaries, as Fortran
304 * compilers (e.g. gfortran) may assume that. Malloc gives suitable
305 * alignment, so we are here careful not to mess it up.
306 */
307 tx = wa;
308 ty = tx + nmax;
309 c = ty + nmax;
310 wrk1 = c + lcest;
311 iwrk = (F_INTint *)(wrk1 + lwrk1);
312 wrk2 = ((double *)iwrk) + kwrk;
313 if (iopt) {
314 ap_tx = (PyArrayObject *)PyArray_ContiguousFromObject(tx_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(tx_py, (*(PyArray_Descr * (*)(
int)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 |
0x0400))) | 0x0040, ((void*)0));
315 ap_ty = (PyArrayObject *)PyArray_ContiguousFromObject(ty_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(ty_py, (*(PyArray_Descr * (*)(
int)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 |
0x0400))) | 0x0040, ((void*)0));
316 if (ap_tx == NULL((void*)0) || ap_ty == NULL((void*)0)) {
317 goto fail;
318 }
319 nx = nxo = PyArray_DIMS(ap_tx)[0];
320 ny = nyo = PyArray_DIMS(ap_ty)[0];
321 memcpy(tx, PyArray_DATA(ap_tx), nx*sizeof(double));
322 memcpy(ty, PyArray_DATA(ap_ty), ny*sizeof(double));
323 }
324 if (iopt==1) {
325 lc = (nx - kx - 1)*(ny - ky - 1);
326 memcpy(wrk1, PyArray_DATA(ap_wrk), lc*sizeof(double));
327 /*memcpy(iwrk,PyArray_DATA(ap_iwrk),n*sizeof(int));*/
328 }
329 SURFITsurfit_(&iopt, &m, x, y, z, w, &xb, &xe, &yb, &ye, &kx, &ky,
330 &s, &nxest, &nyest, &nmax, &eps, &nx, tx, &ny, ty,
331 c, &fp, wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, &ier);
332 i = 0;
333 while ((ier > 10) && (i++ < 5)) {
334 lwrk2 = ier;
335 if ((wrk2 = malloc(lwrk2*sizeof(double))) == NULL((void*)0)) {
336 PyErr_NoMemory();
337 goto fail;
338 }
339 SURFITsurfit_(&iopt, &m, x, y, z, w, &xb, &xe, &yb, &ye, &kx, &ky,
340 &s, &nxest, &nyest, &nmax, &eps, &nx, tx, &ny, ty,
341 c, &fp, wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, &ier);
342 free(wrk2);
343 }
344 if (ier == 10) {
345 PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
346 goto fail;
347 }
348 lc = (nx - kx - 1)*(ny - ky - 1);
349 Py_XDECREF(ap_tx)_Py_XDECREF(((PyObject*)(ap_tx)));
350 Py_XDECREF(ap_ty)_Py_XDECREF(((PyObject*)(ap_ty)));
351 dims[0] = nx;
352 ap_tx = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
353 dims[0] = ny;
354 ap_ty = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
355 dims[0] = lc;
356 ap_c = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
357 if (ap_tx == NULL((void*)0)
358 || ap_ty == NULL((void*)0)
359 || ap_c == NULL((void*)0)) {
360 goto fail;
361 }
362 if ((iopt == 0)||(nx > nxo)||(ny > nyo)) {
363 Py_XDECREF(ap_wrk)_Py_XDECREF(((PyObject*)(ap_wrk)));
364 dims[0] = lc;
365 ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
366 if (ap_wrk == NULL((void*)0)) {
367 goto fail;
368 }
369 /*ap_iwrk = (PyArrayObject *)PyArray_SimpleNew(1,&n,F_INT_NPY);*/
370 }
371 if (PyArray_DIMS(ap_wrk)[0] < lc) {
372 Py_XDECREF(ap_wrk)_Py_XDECREF(((PyObject*)(ap_wrk)));
373 dims[0] = lc;
374 ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
375 if (ap_wrk == NULL((void*)0)) {
376 goto fail;
377 }
378 }
379 memcpy(PyArray_DATA(ap_tx), tx, nx*sizeof(double));
380 memcpy(PyArray_DATA(ap_ty), ty, ny*sizeof(double));
381 memcpy(PyArray_DATA(ap_c), c, lc*sizeof(double));
382 memcpy(PyArray_DATA(ap_wrk), wrk1, lc*sizeof(double));
383 /*memcpy(PyArray_DATA(ap_iwrk),iwrk,n*sizeof(int));*/
384 free(wa);
385 Py_DECREF(ap_x)_Py_DECREF(((PyObject*)(ap_x)));
386 Py_DECREF(ap_y)_Py_DECREF(((PyObject*)(ap_y)));
387 Py_DECREF(ap_z)_Py_DECREF(((PyObject*)(ap_z)));
388 Py_DECREF(ap_w)_Py_DECREF(((PyObject*)(ap_w)));
389 return Py_BuildValue("NNN{s:N,s:i,s:d}",PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_tx),
390 PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_ty),PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_c),
391 "wrk",PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_wrk),
392 "ier",ier,"fp",fp);
393
394fail:
395 free(wa);
396 Py_XDECREF(ap_x)_Py_XDECREF(((PyObject*)(ap_x)));
397 Py_XDECREF(ap_y)_Py_XDECREF(((PyObject*)(ap_y)));
398 Py_XDECREF(ap_z)_Py_XDECREF(((PyObject*)(ap_z)));
399 Py_XDECREF(ap_w)_Py_XDECREF(((PyObject*)(ap_w)));
400 Py_XDECREF(ap_tx)_Py_XDECREF(((PyObject*)(ap_tx)));
401 Py_XDECREF(ap_ty)_Py_XDECREF(((PyObject*)(ap_ty)));
402 Py_XDECREF(ap_wrk)_Py_XDECREF(((PyObject*)(ap_wrk)));
403 /*Py_XDECREF(ap_iwrk);*/
404 if (!PyErr_Occurred()) {
405 PyErr_SetString(PyExc_ValueError,
406 "An error occurred.");
407 }
408 return NULL((void*)0);
409}
410
411
412static char doc_parcur[] = " [t,c,o] = _parcur(x,w,u,ub,ue,k,iopt,ipar,s,t,nest,wrk,iwrk,per)";
413static PyObject *
414fitpack_parcur(PyObject *dummy, PyObject *args)
415{
416 F_INTint k, iopt, ipar, nest, *iwrk, idim, m, mx, no=0, nc, ier, lwa, lwrk, i, per;
417 F_INTint n=0, lc;
418 npy_intp dims[1];
419 double *x, *w, *u, *c, *t, *wrk, *wa=NULL((void*)0), ub, ue, fp, s;
420 PyObject *x_py = NULL((void*)0), *u_py = NULL((void*)0), *w_py = NULL((void*)0), *t_py = NULL((void*)0);
421 PyObject *wrk_py=NULL((void*)0), *iwrk_py=NULL((void*)0);
422 PyArrayObject *ap_x = NULL((void*)0), *ap_u = NULL((void*)0), *ap_w = NULL((void*)0), *ap_t = NULL((void*)0), *ap_c = NULL((void*)0);
423 PyArrayObject *ap_wrk = NULL((void*)0), *ap_iwrk = NULL((void*)0);
424
425 if (!PyArg_ParseTuple(args, ("OOOdd" F_INT_PYFMT"i" F_INT_PYFMT"i" F_INT_PYFMT"i"
426 "dO" F_INT_PYFMT"i" "OO" F_INT_PYFMT"i"),
427 &x_py, &w_py, &u_py, &ub, &ue, &k, &iopt, &ipar,
428 &s, &t_py, &nest, &wrk_py, &iwrk_py, &per)) {
429 return NULL((void*)0);
430 }
431 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
432 ap_u = (PyArrayObject *)PyArray_ContiguousFromObject(u_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(u_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
433 ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(w_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
434 ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(wrk_py, (*(PyArray_Descr * (*)
(int)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100
| 0x0400))) | 0x0040, ((void*)0));
435 ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, F_INT_NPY, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(iwrk_py, (*(PyArray_Descr * (*
)(int)) PyArray_API[45])(NPY_INT), 0, 1, ((0x0001 | (0x0100 |
0x0400))) | 0x0040, ((void*)0));
436 if (ap_x == NULL((void*)0)
437 || ap_u == NULL((void*)0)
438 || ap_w == NULL((void*)0)
439 || ap_wrk == NULL((void*)0)
440 || ap_iwrk == NULL((void*)0)) {
441 goto fail;
442 }
443 x = (double *) PyArray_DATA(ap_x);
444 u = (double *) PyArray_DATA(ap_u);
445 w = (double *) PyArray_DATA(ap_w);
446 m = PyArray_DIMS(ap_w)[0];
447 mx = PyArray_DIMS(ap_x)[0];
448 idim = mx/m;
449 if (per) {
450 lwrk = m*(k + 1) + nest*(7 + idim + 5*k);
451 }
452 else {
453 lwrk = m*(k + 1) + nest*(6 + idim + 3*k);
454 }
455 nc = idim*nest;
456 lwa = nc + 2*nest + lwrk;
457 if ((wa = malloc(lwa*sizeof(double))) == NULL((void*)0)) {
458 PyErr_NoMemory();
459 goto fail;
460 }
461 t = wa;
462 c = t + nest;
463 wrk = c + nc;
464 iwrk = (F_INTint *)(wrk + lwrk);
465 if (iopt) {
466 ap_t=(PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(t_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
467 if (ap_t == NULL((void*)0)) {
468 goto fail;
469 }
470 n = no = PyArray_DIMS(ap_t)[0];
471 memcpy(t, PyArray_DATA(ap_t), n*sizeof(double));
472 Py_DECREF(ap_t)_Py_DECREF(((PyObject*)(ap_t)));
473 ap_t = NULL((void*)0);
474 }
475 if (iopt == 1) {
476 memcpy(wrk, PyArray_DATA(ap_wrk), n*sizeof(double));
477 memcpy(iwrk, PyArray_DATA(ap_iwrk), n*sizeof(F_INTint));
478 }
479 if (per) {
480 CLOCURclocur_(&iopt, &ipar, &idim, &m, u, &mx, x, w, &k, &s, &nest,
481 &n, t, &nc, c, &fp, wrk, &lwrk, iwrk, &ier);
482 }
483 else {
484 PARCURparcur_(&iopt, &ipar, &idim, &m, u, &mx, x, w, &ub, &ue, &k,
485 &s, &nest, &n, t, &nc, c, &fp, wrk, &lwrk, iwrk, &ier);
486 }
487 if (ier == 10) {
488 PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
489 goto fail;
490 }
491 if (ier > 0 && n == 0) {
492 n = 1;
493 }
494 lc = (n - k - 1)*idim;
495 dims[0] = n;
496 ap_t = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
497 dims[0] = lc;
498 ap_c = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
499 if (ap_t == NULL((void*)0) || ap_c == NULL((void*)0)) {
500 goto fail;
501 }
502 if (iopt != 1|| n > no) {
503 Py_XDECREF(ap_wrk)_Py_XDECREF(((PyObject*)(ap_wrk)));
504 ap_wrk = NULL((void*)0);
505 Py_XDECREF(ap_iwrk)_Py_XDECREF(((PyObject*)(ap_iwrk)));
506 ap_iwrk = NULL((void*)0);
507
508 dims[0] = n;
509 ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
510 if (ap_wrk == NULL((void*)0)) {
511 goto fail;
512 }
513 ap_iwrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, F_INT_NPY)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_INT
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
514 if (ap_iwrk == NULL((void*)0)) {
515 goto fail;
516 }
517 }
518 memcpy(PyArray_DATA(ap_t), t, n*sizeof(double));
519 for (i = 0; i < idim; i++)
520 memcpy((double *)PyArray_DATA(ap_c) + i*(n - k - 1), c + i*n, (n - k - 1)*sizeof(double));
521 memcpy(PyArray_DATA(ap_wrk), wrk, n*sizeof(double));
522 memcpy(PyArray_DATA(ap_iwrk), iwrk, n*sizeof(F_INTint));
523 free(wa);
524 Py_DECREF(ap_x)_Py_DECREF(((PyObject*)(ap_x)));
525 Py_DECREF(ap_w)_Py_DECREF(((PyObject*)(ap_w)));
526 return Py_BuildValue(("NN{s:N,s:d,s:d,s:N,s:N,s:" F_INT_PYFMT"i" ",s:d}"), PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_t),
527 PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_c), "u", PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_u), "ub", ub, "ue", ue,
528 "wrk", PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_wrk), "iwrk", PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_iwrk),
529 "ier", ier, "fp",fp);
530fail:
531 free(wa);
532 Py_XDECREF(ap_x)_Py_XDECREF(((PyObject*)(ap_x)));
533 Py_XDECREF(ap_u)_Py_XDECREF(((PyObject*)(ap_u)));
534 Py_XDECREF(ap_w)_Py_XDECREF(((PyObject*)(ap_w)));
535 Py_XDECREF(ap_t)_Py_XDECREF(((PyObject*)(ap_t)));
536 Py_XDECREF(ap_wrk)_Py_XDECREF(((PyObject*)(ap_wrk)));
537 Py_XDECREF(ap_iwrk)_Py_XDECREF(((PyObject*)(ap_iwrk)));
538 return NULL((void*)0);
539}
540
541static char doc_curfit[] = " [t,c,o] = _curfit(x,y,w,xb,xe,k,iopt,s,t,nest,wrk,iwrk,per)";
542static PyObject *
543fitpack_curfit(PyObject *dummy, PyObject *args)
544{
545 F_INTint iopt, m, k, nest, lwrk, *iwrk, ier, lwa, no=0, per;
546 F_INTint n, lc;
547 npy_intp dims[1];
548 double *x, *y, *w, xb, xe, s, *t, *c, fp, *wrk, *wa = NULL((void*)0);
549 PyArrayObject *ap_x = NULL((void*)0), *ap_y = NULL((void*)0), *ap_w = NULL((void*)0), *ap_t = NULL((void*)0), *ap_c = NULL((void*)0);
550 PyArrayObject *ap_wrk = NULL((void*)0), *ap_iwrk = NULL((void*)0);
551 PyObject *x_py = NULL((void*)0), *y_py = NULL((void*)0), *w_py = NULL((void*)0), *t_py = NULL((void*)0);
552 PyObject *wrk_py=NULL((void*)0), *iwrk_py=NULL((void*)0);
553
554 if (!PyArg_ParseTuple(args, ("OOOdd" F_INT_PYFMT"i" F_INT_PYFMT"i"
555 "dO" F_INT_PYFMT"i" "OO" F_INT_PYFMT"i"),
556 &x_py, &y_py, &w_py, &xb, &xe, &k, &iopt,
557 &s, &t_py, &nest, &wrk_py, &iwrk_py, &per)) {
558 return NULL((void*)0);
559 }
560 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
561 ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(y_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
562 ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(w_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
563 ap_wrk = (PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(wrk_py, (*(PyArray_Descr * (*)
(int)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100
| 0x0400))) | 0x0040, ((void*)0));
564 ap_iwrk = (PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, F_INT_NPY, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(iwrk_py, (*(PyArray_Descr * (*
)(int)) PyArray_API[45])(NPY_INT), 0, 1, ((0x0001 | (0x0100 |
0x0400))) | 0x0040, ((void*)0));
565 if (ap_x == NULL((void*)0)
566 || ap_y == NULL((void*)0)
567 || ap_w == NULL((void*)0)
568 || ap_wrk == NULL((void*)0)
569 || ap_iwrk == NULL((void*)0)) {
570 goto fail;
571 }
572 x = (double *) PyArray_DATA(ap_x);
573 y = (double *) PyArray_DATA(ap_y);
574 w = (double *) PyArray_DATA(ap_w);
575 m = PyArray_DIMS(ap_x)[0];
576 if (per) {
577 lwrk = m*(k + 1) + nest*(8 + 5*k);
578 }
579 else {
580 lwrk = m*(k + 1) + nest*(7 + 3*k);
581 }
582 lwa = 3*nest + lwrk;
583 if ((wa = malloc(lwa*sizeof(double))) == NULL((void*)0)) {
584 PyErr_NoMemory();
585 goto fail;
586 }
587 t = wa;
588 c = t + nest;
589 wrk = c + nest;
590 iwrk = (F_INTint *)(wrk + lwrk);
591 if (iopt) {
592 ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(t_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
593 if (ap_t == NULL((void*)0)) {
594 goto fail;
595 }
596 n = no = PyArray_DIMS(ap_t)[0];
597 memcpy(t, PyArray_DATA(ap_t), n*sizeof(double));
598 }
599 if (iopt == 1) {
600 memcpy(wrk, PyArray_DATA(ap_wrk), n*sizeof(double));
601 memcpy(iwrk, PyArray_DATA(ap_iwrk), n*sizeof(F_INTint));
602 }
603 if (per)
604 PERCURpercur_(&iopt, &m, x, y, w, &k, &s, &nest, &n, t, c, &fp, wrk,
605 &lwrk, iwrk, &ier);
606 else {
607 CURFITcurfit_(&iopt, &m, x, y, w, &xb, &xe, &k, &s, &nest, &n, t, c,
608 &fp, wrk, &lwrk, iwrk, &ier);
609 }
610 if (ier == 10) {
611 PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
612 goto fail;
613 }
614 lc = n - k - 1;
615 if (!iopt) {
616 dims[0] = n;
617 ap_t = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
618 if (ap_t == NULL((void*)0)) {
619 goto fail;
620 }
621 }
622 dims[0] = lc;
623 ap_c = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
624 if (ap_c == NULL((void*)0)) {
625 goto fail;
626 }
627 if (iopt == 0 || n > no) {
628 Py_XDECREF(ap_wrk)_Py_XDECREF(((PyObject*)(ap_wrk)));
629 Py_XDECREF(ap_iwrk)_Py_XDECREF(((PyObject*)(ap_iwrk)));
630 dims[0] = n;
631 ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
632 ap_iwrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, F_INT_NPY)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_INT
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
633 if (ap_wrk == NULL((void*)0) || ap_iwrk == NULL((void*)0)) {
634 goto fail;
635 }
636 }
637 memcpy(PyArray_DATA(ap_t), t, n*sizeof(double));
638 memcpy(PyArray_DATA(ap_c), c, lc*sizeof(double));
639 memcpy(PyArray_DATA(ap_wrk), wrk, n*sizeof(double));
640 memcpy(PyArray_DATA(ap_iwrk), iwrk, n*sizeof(F_INTint));
641 free(wa);
642 Py_DECREF(ap_x)_Py_DECREF(((PyObject*)(ap_x)));
643 Py_DECREF(ap_y)_Py_DECREF(((PyObject*)(ap_y)));
644 Py_DECREF(ap_w)_Py_DECREF(((PyObject*)(ap_w)));
645 return Py_BuildValue(("NN{s:N,s:N,s:" F_INT_PYFMT"i" ",s:d}"), PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_t),
646 PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_c), "wrk", PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_wrk),
647 "iwrk", PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_iwrk), "ier", ier, "fp", fp);
648
649fail:
650 free(wa);
651 Py_XDECREF(ap_x)_Py_XDECREF(((PyObject*)(ap_x)));
652 Py_XDECREF(ap_y)_Py_XDECREF(((PyObject*)(ap_y)));
653 Py_XDECREF(ap_w)_Py_XDECREF(((PyObject*)(ap_w)));
654 Py_XDECREF(ap_t)_Py_XDECREF(((PyObject*)(ap_t)));
655 Py_XDECREF(ap_wrk)_Py_XDECREF(((PyObject*)(ap_wrk)));
656 Py_XDECREF(ap_iwrk)_Py_XDECREF(((PyObject*)(ap_iwrk)));
657 return NULL((void*)0);
658}
659
660static char doc_spl_[] = " [y,ier] = _spl_(x,nu,t,c,k,e)";
661static PyObject *
662fitpack_spl_(PyObject *dummy, PyObject *args)
663{
664 F_INTint n, nu, ier, k, m, e=0;
665 npy_intp dims[1];
666 double *x, *y, *t, *c, *wrk = NULL((void*)0);
667 PyArrayObject *ap_x = NULL((void*)0), *ap_y = NULL((void*)0), *ap_t = NULL((void*)0), *ap_c = NULL((void*)0);
668 PyObject *x_py = NULL((void*)0), *t_py = NULL((void*)0), *c_py = NULL((void*)0);
669
670 if (!PyArg_ParseTuple(args, ("O" F_INT_PYFMT"i" "OO" F_INT_PYFMT"i" F_INT_PYFMT"i"),
671 &x_py, &nu, &t_py, &c_py, &k, &e)) {
672 return NULL((void*)0);
673 }
674 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
675 ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(t_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
676 ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(c_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
677 if ((ap_x == NULL((void*)0) || ap_t == NULL((void*)0) || ap_c == NULL((void*)0))) {
678 goto fail;
679 }
680 x = (double *)PyArray_DATA(ap_x);
681 m = PyArray_DIMS(ap_x)[0];
682 t = (double *)PyArray_DATA(ap_t);
683 c = (double *)PyArray_DATA(ap_c);
684 n = PyArray_DIMS(ap_t)[0];
685 dims[0] = m;
686 ap_y = (PyArrayObject *)PyArray_SimpleNew(1,dims,NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
687 if (ap_y == NULL((void*)0)) {
688 goto fail;
689 }
690 y = (double *)PyArray_DATA(ap_y);
691 if ((wrk = malloc(n*sizeof(double))) == NULL((void*)0)) {
692 PyErr_NoMemory();
693 goto fail;
694 }
695 if (nu) {
696 SPLDERsplder_(t, &n, c, &k, &nu, x, y, &m, &e, wrk, &ier);
697 }
698 else {
699 SPLEVsplev_(t, &n, c, &k, x, y, &m, &e, &ier);
700 }
701 free(wrk);
702 Py_DECREF(ap_x)_Py_DECREF(((PyObject*)(ap_x)));
703 Py_DECREF(ap_c)_Py_DECREF(((PyObject*)(ap_c)));
704 Py_DECREF(ap_t)_Py_DECREF(((PyObject*)(ap_t)));
705 return Py_BuildValue(("N" F_INT_PYFMT"i"), PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_y), ier);
706
707fail:
708 free(wrk);
709 Py_XDECREF(ap_x)_Py_XDECREF(((PyObject*)(ap_x)));
710 Py_XDECREF(ap_c)_Py_XDECREF(((PyObject*)(ap_c)));
711 Py_XDECREF(ap_t)_Py_XDECREF(((PyObject*)(ap_t)));
712 return NULL((void*)0);
713}
714
715static char doc_splint[] = " [aint,wrk] = _splint(t,c,k,a,b)";
716static PyObject *
717fitpack_splint(PyObject *dummy, PyObject *args)
718{
719 F_INTint k, n;
720 npy_intp dims[1];
721 double *t, *c, *wrk = NULL((void*)0), a, b, aint;
722 PyArrayObject *ap_t = NULL((void*)0), *ap_c = NULL((void*)0);
723 PyArrayObject *ap_wrk = NULL((void*)0);
724 PyObject *t_py = NULL((void*)0), *c_py = NULL((void*)0);
725
726 if (!PyArg_ParseTuple(args, ("OO" F_INT_PYFMT"i" "dd"),&t_py,&c_py,&k,&a,&b)) {
727 return NULL((void*)0);
728 }
729 ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(t_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
730 ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(c_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
731 if ((ap_t == NULL((void*)0) || ap_c == NULL((void*)0))) {
732 goto fail;
733 }
734 t = (double *)PyArray_DATA(ap_t);
735 c = (double *)PyArray_DATA(ap_c);
736 n = PyArray_DIMS(ap_t)[0];
737 dims[0] = n;
738 ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
739 if (ap_wrk == NULL((void*)0)) {
740 goto fail;
741 }
742 wrk = (double *)PyArray_DATA(ap_wrk);
743 aint = SPLINTsplint_(t,&n,c,&k,&a,&b,wrk);
744 Py_DECREF(ap_c)_Py_DECREF(((PyObject*)(ap_c)));
745 Py_DECREF(ap_t)_Py_DECREF(((PyObject*)(ap_t)));
746 return Py_BuildValue("dN", aint, PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_wrk));
747
748fail:
749 Py_XDECREF(ap_c)_Py_XDECREF(((PyObject*)(ap_c)));
750 Py_XDECREF(ap_t)_Py_XDECREF(((PyObject*)(ap_t)));
751 return NULL((void*)0);
752}
753
754static char doc_sproot[] = " [z,ier] = _sproot(t,c,k,mest)";
755static PyObject *
756fitpack_sproot(PyObject *dummy, PyObject *args)
757{
758 F_INTint n, k, m, mest, ier;
759 npy_intp dims[1];
760 double *t, *c, *z = NULL((void*)0);
761 PyArrayObject *ap_t = NULL((void*)0), *ap_c = NULL((void*)0);
762 PyArrayObject *ap_z = NULL((void*)0);
763 PyObject *t_py = NULL((void*)0), *c_py = NULL((void*)0);
764
765 if (!PyArg_ParseTuple(args, ("OO" F_INT_PYFMT"i" F_INT_PYFMT"i"),
766 &t_py,&c_py,&k,&mest)) {
767 return NULL((void*)0);
768 }
769 ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(t_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
770 ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(c_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
771 if ((ap_t == NULL((void*)0) || ap_c == NULL((void*)0))) {
772 goto fail;
773 }
774 t = (double *)PyArray_DATA(ap_t);
775 c = (double *)PyArray_DATA(ap_c);
776 n = PyArray_DIMS(ap_t)[0];
777 if ((z = malloc(mest*sizeof(double))) == NULL((void*)0)) {
778 PyErr_NoMemory();
779 goto fail;
780 }
781 m = 0;
782 SPROOTsproot_(t,&n,c,z,&mest,&m,&ier);
783 if (ier==10) {
784 m = 0;
785 }
786 dims[0] = m;
787 ap_z = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
788 if (ap_z == NULL((void*)0)) {
789 goto fail;
790 }
791 memcpy(PyArray_DATA(ap_z), z, m*sizeof(double));
792 free(z);
793 Py_DECREF(ap_c)_Py_DECREF(((PyObject*)(ap_c)));
794 Py_DECREF(ap_t)_Py_DECREF(((PyObject*)(ap_t)));
795 return Py_BuildValue(("N" F_INT_PYFMT"i"), PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_z), ier);
796
797fail:
798 free(z);
799 Py_XDECREF(ap_c)_Py_XDECREF(((PyObject*)(ap_c)));
800 Py_XDECREF(ap_t)_Py_XDECREF(((PyObject*)(ap_t)));
801 return NULL((void*)0);
802}
803
804static char doc_spalde[] = " [d,ier] = _spalde(t,c,k,x)";
805static PyObject *
806fitpack_spalde(PyObject *dummy, PyObject *args)
807{
808 F_INTint n, k, ier, k1;
809 npy_intp dims[1];
810 double *t, *c, *d = NULL((void*)0), x;
811 PyArrayObject *ap_t = NULL((void*)0), *ap_c = NULL((void*)0), *ap_d = NULL((void*)0);
812 PyObject *t_py = NULL((void*)0), *c_py = NULL((void*)0);
813
814 if (!PyArg_ParseTuple(args, ("OO" F_INT_PYFMT"i" "d"),
815 &t_py,&c_py,&k,&x)) {
816 return NULL((void*)0);
817 }
818 ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(t_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
819 ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(c_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
820 if ((ap_t == NULL((void*)0) || ap_c == NULL((void*)0))) {
821 goto fail;
822 }
823 t = (double *)PyArray_DATA(ap_t);
824 c = (double *)PyArray_DATA(ap_c);
825 n = PyArray_DIMS(ap_t)[0];
826 k1 = k + 1;
827 dims[0] = k1;
828 ap_d = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
829 if (ap_d == NULL((void*)0)) {
830 goto fail;
831 }
832 d = (double *)PyArray_DATA(ap_d);
833 SPALDEspalde_(t, &n, c, &k1, &x, d, &ier);
834 Py_DECREF(ap_c)_Py_DECREF(((PyObject*)(ap_c)));
835 Py_DECREF(ap_t)_Py_DECREF(((PyObject*)(ap_t)));
836 return Py_BuildValue(("N" F_INT_PYFMT"i"), PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_d), ier);
837
838fail:
839 Py_XDECREF(ap_c)_Py_XDECREF(((PyObject*)(ap_c)));
840 Py_XDECREF(ap_t)_Py_XDECREF(((PyObject*)(ap_t)));
841 return NULL((void*)0);
842}
843
844static char doc_insert[] = " [tt,cc,ier] = _insert(iopt,t,c,k,x,m)";
845static PyObject *
846fitpack_insert(PyObject *dummy, PyObject*args)
847{
848 F_INTint iopt, n, nn, k, ier, m, nest;
849 npy_intp dims[1];
850 double x;
851 double *t_in, *c_in, *t_out, *c_out, *t_buf = NULL((void*)0), *c_buf = NULL((void*)0), *p;
852 double *t1, *t2, *c1, *c2;
853 PyArrayObject *ap_t_in = NULL((void*)0), *ap_c_in = NULL((void*)0), *ap_t_out = NULL((void*)0), *ap_c_out = NULL((void*)0);
854 PyObject *t_py = NULL((void*)0), *c_py = NULL((void*)0);
855 PyObject *ret = NULL((void*)0);
856
857 if (!PyArg_ParseTuple(args,(F_INT_PYFMT"i" "OO" F_INT_PYFMT"i" "d" F_INT_PYFMT"i"),
858 &iopt,&t_py,&c_py,&k, &x, &m)) {
859 return NULL((void*)0);
860 }
861 ap_t_in = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(t_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
862 ap_c_in = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(c_py, (*(PyArray_Descr * (*)(int
)) PyArray_API[45])(NPY_DOUBLE), 0, 1, ((0x0001 | (0x0100 | 0x0400
))) | 0x0040, ((void*)0));
863 if (ap_t_in == NULL((void*)0) || ap_c_in == NULL((void*)0)) {
864 goto fail;
865 }
866 t_in = (double *)PyArray_DATA(ap_t_in);
867 c_in = (double *)PyArray_DATA(ap_c_in);
868 n = PyArray_DIMS(ap_t_in)[0];
869 nest = n + m;
870 dims[0] = nest;
871 ap_t_out = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
872 ap_c_out = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE)(*(PyObject * (*)(PyTypeObject *, int, npy_intp const *, int,
npy_intp const *, void *, int, int, PyObject *)) PyArray_API
[93])(&(*(PyTypeObject *)PyArray_API[2]), 1, dims, NPY_DOUBLE
, ((void*)0), ((void*)0), 0, 0, ((void*)0));
873 if (ap_t_out == NULL((void*)0) || ap_c_out == NULL((void*)0)) {
874 goto fail;
875 }
876 t_out = (double *)PyArray_DATA(ap_t_out);
877 c_out = (double *)PyArray_DATA(ap_c_out);
878
879 /*
880 * Call the INSERT routine m times to insert m-multiplicity knot, ie.:
881 *
882 * for _ in range(n, nest):
883 * t, c = INSERT(t, c)
884 * return t, c
885 *
886 * We need to ensure that input and output buffers given to INSERT routine
887 * do not point to same memory, which is not allowed by Fortran. For this,
888 * we use temporary storage, and cycle between it and the output buffers.
889 */
890 t2 = t_in;
891 c2 = c_in;
892 t1 = t_out;
893 c1 = c_out;
894
895 for ( ; n < nest; n++) {
896 /* Swap buffers */
897 p = t2; t2 = t1; t1 = p;
898 p = c2; c2 = c1; c1 = p;
899
900 /* Allocate temporary buffer (needed for m > 1) */
901 if (t2 == t_in) {
902 if (t_buf == NULL((void*)0)) {
903 t_buf = calloc(nest, sizeof(double));
904 c_buf = calloc(nest, sizeof(double));
905 if (t_buf == NULL((void*)0) || c_buf == NULL((void*)0)) {
906 PyErr_NoMemory();
907 goto fail;
908 }
909 }
910 t2 = t_buf;
911 c2 = c_buf;
912 }
913
914 INSERTinsert_(&iopt, t1, &n, c1, &k, &x, t2, &nn, c2, &nest, &ier);
915
916 if (ier) {
917 break;
918 }
919 }
920
921 /* Ensure output ends up in output buffers */
922 if (t2 != t_out) {
923 memcpy(t_out, t2, nest * sizeof(double));
924 memcpy(c_out, c2, nest * sizeof(double));
925 }
926
927 Py_DECREF(ap_c_in)_Py_DECREF(((PyObject*)(ap_c_in)));
928 Py_DECREF(ap_t_in)_Py_DECREF(((PyObject*)(ap_t_in)));
929 free(t_buf);
930 free(c_buf);
931 ret = Py_BuildValue(("NN" F_INT_PYFMT"i"), PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_t_out), PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(ap_c_out), ier);
932 return ret;
933
934fail:
935 Py_XDECREF(ap_c_out)_Py_XDECREF(((PyObject*)(ap_c_out)));
936 Py_XDECREF(ap_t_out)_Py_XDECREF(((PyObject*)(ap_t_out)));
937 Py_XDECREF(ap_c_in)_Py_XDECREF(((PyObject*)(ap_c_in)));
938 Py_XDECREF(ap_t_in)_Py_XDECREF(((PyObject*)(ap_t_in)));
939 free(t_buf);
940 free(c_buf);
941 return NULL((void*)0);
942}
943
944
945/*
946 * Given a set of (N+1) samples: A default set of knots is constructed
947 * using the samples xk plus 2*(K-1) additional knots where
948 * K = max(order,1) and the knots are chosen so that distances
949 * are symmetric around the first and last samples: x_0 and x_N.
950 *
951 * There should be a vector of N+K coefficients for the spline
952 * curve in coef. These coefficients form the curve as
953 *
954 * s(x) = sum(c_j B_{j,K}(x), j=-K..N-1)
955 *
956 * The spline function is evaluated at all points xx.
957 * The approximation interval is from xk[0] to xk[-1]
958 * Any xx outside that interval is set automatically to 0.0
959 */
960static char doc_bspleval[] = "y = _bspleval(xx,xk,coef,k,{deriv (0)})\n"
961"\n"
962"The spline is defined by the approximation interval xk[0] to xk[-1],\n"
963"the length of xk (N+1), the order of the spline, k, and \n"
964"the number of coeficients N+k. The coefficients range from xk_{-K}\n"
965"to xk_{N-1} inclusive and are all the coefficients needed to define\n"
966"an arbitrary spline of order k, on the given approximation interval\n"
967"\n"
968"Extra knot points are internally added using knot-point symmetry \n"
969"around xk[0] and xk[-1]";
970
971static PyObject *
972_bspleval(PyObject *dummy, PyObject *args)
973{
974 int k, kk, N, i, ell, dk, deriv = 0;
975 PyObject *xx_py = NULL((void*)0), *coef_py = NULL((void*)0), *x_i_py = NULL((void*)0);
976 PyArrayObject *xx = NULL((void*)0), *coef = NULL((void*)0), *x_i = NULL((void*)0), *yy = NULL((void*)0);
977 PyArrayIterObject *xx_iter;
978 double *t = NULL((void*)0), *h = NULL((void*)0), *ptr;
979 double x0, xN, xN1, arg, sp, cval;
980
981 if (!PyArg_ParseTuple(args, "OOOi|i", &xx_py, &x_i_py, &coef_py, &k, &deriv)) {
982 return NULL((void*)0);
983 }
984 if (k < 0) {
985 PyErr_Format(PyExc_ValueError,
986 "order (%d) must be >=0", k);
987 return NULL((void*)0);
988 }
989 if (deriv > k) {
990 PyErr_Format(PyExc_ValueError,
991 "derivative (%d) must be <= order (%d)", deriv, k);
992 return NULL((void*)0);
993 }
994 kk = k;
995 if (k == 0) {
996 kk = 1;
997 }
998 dk = (k == 0 ? 0 : 1);
999 x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ARRAY_ALIGNED)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_i_py, (*(PyArray_Descr * (*)
(int)) PyArray_API[45])(NPY_DOUBLE), 1, 1, (((0x0100) & 0x0020
) ? (0x0100) | ((0x0001 | (0x0100 | 0x0400))) : (0x0100)), ((
void*)0));
1000 coef = (PyArrayObject *)PyArray_FROMANY(coef_py, NPY_DOUBLE, 1, 1, NPY_ARRAY_ALIGNED)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(coef_py, (*(PyArray_Descr * (*
)(int)) PyArray_API[45])(NPY_DOUBLE), 1, 1, (((0x0100) & 0x0020
) ? (0x0100) | ((0x0001 | (0x0100 | 0x0400))) : (0x0100)), ((
void*)0));
1001 xx = (PyArrayObject *)PyArray_FROMANY(xx_py, NPY_DOUBLE, 0, 0, NPY_ARRAY_ALIGNED)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(xx_py, (*(PyArray_Descr * (*)(
int)) PyArray_API[45])(NPY_DOUBLE), 0, 0, (((0x0100) & 0x0020
) ? (0x0100) | ((0x0001 | (0x0100 | 0x0400))) : (0x0100)), ((
void*)0));
1002 if (x_i == NULL((void*)0) || coef == NULL((void*)0) || xx == NULL((void*)0)) {
1003 goto fail;
1004 }
1005
1006 N = PyArray_DIM(x_i, 0) - 1;
1007 if (PyArray_DIM(coef, 0) < N + k) {
1008 PyErr_Format(PyExc_ValueError,
1009 "too few coefficients (have %d need at least %d)",
1010 (int)PyArray_DIM(coef, 0), N + k);
1011 goto fail;
1012 }
1013
1014 /* create output values */
1015 yy = (PyArrayObject *)PyArray_EMPTY(PyArray_NDIM(xx), PyArray_DIMS(xx), NPY_DOUBLE, 0)(*(PyObject * (*)(int, npy_intp const *, PyArray_Descr *, int
)) PyArray_API[184])(PyArray_NDIM(xx), PyArray_DIMS(xx), (*(PyArray_Descr
* (*)(int)) PyArray_API[45])(NPY_DOUBLE), 0);
1016 if (yy == NULL((void*)0)) {
1017 goto fail;
1018 }
1019 /*
1020 * create dummy knot array with new knots inserted at the end
1021 * selected as mirror symmetric versions of the old knots
1022 */
1023 t = malloc(sizeof(double)*(N + 2*kk - 1));
1024 if (t == NULL((void*)0)) {
1025 PyErr_NoMemory();
1026 goto fail;
1027 }
1028 x0 = *((double *)PyArray_DATA(x_i));
1029 xN = *((double *)PyArray_DATA(x_i) + N);
1030 for (i = 0; i < kk - 1; i++) { /* fill in ends if kk > 1*/
1031 t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i, kk - 1 - i)((void *)(PyArray_BYTES(x_i) + (kk - 1 - i)*PyArray_STRIDES(x_i
)[0]))));
1032 t[kk+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i, N - 1 - i)((void *)(PyArray_BYTES(x_i) + (N - 1 - i)*PyArray_STRIDES(x_i
)[0]))));
1033 }
1034 ptr = t + (kk - 1);
1035 for (i = 0; i <= N; i++) {
1036 *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)((void *)(PyArray_BYTES(x_i) + (i)*PyArray_STRIDES(x_i)[0]))));
1037 }
1038
1039 /*
1040 * Create work array to hold computed non-zero values for
1041 * the spline for a value of x.
1042 */
1043 h = malloc(sizeof(double)*(2*kk+1));
1044 if (h == NULL((void*)0)) {
1045 PyErr_NoMemory();
1046 goto fail;
1047 }
1048
1049 /* Determine the spline for each value of x */
1050 xx_iter = (PyArrayIterObject *)PyArray_IterNew(*(PyObject * (*)(PyObject *)) PyArray_API[98])((PyObject *)xx);
1051 if (xx_iter == NULL((void*)0)) {
1052 goto fail;
1053 }
1054 ptr = PyArray_DATA(yy);
1055
1056 while(PyArray_ITER_NOTDONE(xx_iter)(((PyArrayIterObject *)(xx_iter))->index < ((PyArrayIterObject
*)(xx_iter))->size)) {
1057 arg = *((double *)PyArray_ITER_DATA(xx_iter)((void *)(((PyArrayIterObject *)(xx_iter))->dataptr)));
1058 if ((arg < x0) || (arg > xN)) {
1059 /*
1060 * If we are outside the interpolation region,
1061 * fill with zeros
1062 */
1063 *ptr++ = 0.0;
1064 }
1065 else {
1066 /*
1067 * Find the interval that arg lies between in the set of knots
1068 * t[ell] <= arg < t[ell+1] (last-knot use the previous interval)
1069 */
1070 xN1 = *((double *)PyArray_DATA(x_i) + N-1);
1071 if (arg >= xN1) {
1072 ell = N + kk - 2;
1073 }
1074 else {
1075 ell = kk - 1;
1076 while ((arg > t[ell])) {
1077 ell++;
1078 }
1079 if (arg != t[ell]) {
1080 ell--;
1081 }
1082 }
1083
1084 _deBoor_D(t, arg, k, ell, deriv, h);
1085
1086 sp = 0.0;
1087 for (i = 0; i <= k; i++) {
1088 cval = *((double *)(PyArray_GETPTR1(coef, ell - i + dk)((void *)(PyArray_BYTES(coef) + (ell - i + dk)*PyArray_STRIDES
(coef)[0]))));
1089 sp += cval * h[k - i];
1090 }
1091 *ptr++ = sp;
1092 }
1093 PyArray_ITER_NEXT(xx_iter)do { ((PyArrayIterObject *)(xx_iter))->index++; if (((PyArrayIterObject
*)(xx_iter))->nd_m1 == 0) { do { (((PyArrayIterObject *)(
xx_iter)))->dataptr += ((PyArrayIterObject *)(((PyArrayIterObject
*)(xx_iter))))->strides[0]; (((PyArrayIterObject *)(xx_iter
)))->coordinates[0]++; } while (0); } else if (((PyArrayIterObject
*)(xx_iter))->contiguous) ((PyArrayIterObject *)(xx_iter)
)->dataptr += PyArray_DESCR(((PyArrayIterObject *)(xx_iter
))->ao)->elsize; else if (((PyArrayIterObject *)(xx_iter
))->nd_m1 == 1) { do { if ((((PyArrayIterObject *)(xx_iter
)))->coordinates[1] < (((PyArrayIterObject *)(xx_iter))
)->dims_m1[1]) { (((PyArrayIterObject *)(xx_iter)))->coordinates
[1]++; (((PyArrayIterObject *)(xx_iter)))->dataptr += (((PyArrayIterObject
*)(xx_iter)))->strides[1]; } else { (((PyArrayIterObject *
)(xx_iter)))->coordinates[1] = 0; (((PyArrayIterObject *)(
xx_iter)))->coordinates[0]++; (((PyArrayIterObject *)(xx_iter
)))->dataptr += (((PyArrayIterObject *)(xx_iter)))->strides
[0] - (((PyArrayIterObject *)(xx_iter)))->backstrides[1]; }
} while (0); } else { int __npy_i; for (__npy_i=((PyArrayIterObject
*)(xx_iter))->nd_m1; __npy_i >= 0; __npy_i--) { if (((
PyArrayIterObject *)(xx_iter))->coordinates[__npy_i] < (
(PyArrayIterObject *)(xx_iter))->dims_m1[__npy_i]) { ((PyArrayIterObject
*)(xx_iter))->coordinates[__npy_i]++; ((PyArrayIterObject
*)(xx_iter))->dataptr += ((PyArrayIterObject *)(xx_iter))
->strides[__npy_i]; break; } else { ((PyArrayIterObject *)
(xx_iter))->coordinates[__npy_i] = 0; ((PyArrayIterObject *
)(xx_iter))->dataptr -= ((PyArrayIterObject *)(xx_iter))->
backstrides[__npy_i]; } } } } while (0);
1094 }
1095 Py_DECREF(xx_iter)_Py_DECREF(((PyObject*)(xx_iter)));
1096 Py_DECREF(x_i)_Py_DECREF(((PyObject*)(x_i)));
1097 Py_DECREF(coef)_Py_DECREF(((PyObject*)(coef)));
1098 Py_DECREF(xx)_Py_DECREF(((PyObject*)(xx)));
1099 free(t);
1100 free(h);
1101 return PyArray_Return(*(PyObject * (*)(PyArrayObject *)) PyArray_API[76])(yy);
1102
1103fail:
1104 Py_XDECREF(xx)_Py_XDECREF(((PyObject*)(xx)));
1105 Py_XDECREF(coef)_Py_XDECREF(((PyObject*)(coef)));
1106 Py_XDECREF(x_i)_Py_XDECREF(((PyObject*)(x_i)));
1107 Py_XDECREF(yy)_Py_XDECREF(((PyObject*)(yy)));
1108 free(t);
1109 free(h);
1110 return NULL((void*)0);
1111}
1112
1113
1114/*
1115 * Given a set of (N+1) sample positions:
1116 * Construct the diagonals of the (N+1) x (N+K) matrix that is needed to find
1117 * the coefficients of a spline fit of order K.
1118 * Note that K>=2 because for K=0,1, the coefficients are just the
1119 * sample values themselves.
1120 *
1121 * The equation that expresses the constraints is
1122 *
1123 * s(x_i) = sum(c_j B_{j,K}(x_i), j=-K..N-1) = w_i for i=0..N
1124 *
1125 * This is equivalent to
1126 *
1127 * w = B*c where c.T = [c_{-K}, c{-K+1}, ..., c_{N-1}] and
1128 *
1129 * Therefore B is an (N+1) times (N+K) matrix with entries
1130 *
1131 * B_{j,K}(x_i) for column j=-K..N-1
1132 * and row i=0..N
1133 *
1134 * This routine takes the N+1 sample positions and the order k and
1135 * constructs the banded constraint matrix B (with k non-zero diagonals)
1136 *
1137 * The returned array is (N+1) times (N+K) ready to be either used
1138 * to compute a minimally Kth-order derivative discontinuous spline
1139 * or to be expanded with an additional K-1 constraints to be used in
1140 * an exact spline specification.
1141 */
1142static char doc_bsplmat[] = "B = _bsplmat(order,xk)\n"
1143"Construct the constraint matrix for spline fitting of order k\n"
1144"given sample positions in xk.\n"
1145"\n"
1146"If xk is an integer (N+1), then the result is equivalent to\n"
1147"xk=arange(N+1)+x0 for any value of x0. This produces the\n"
1148"integer-spaced, or cardinal spline matrix a bit faster.";
1149static PyObject *
1150_bsplmat(PyObject *dummy, PyObject *args) {
1151 int k, N, i, numbytes, j, equal;
1152 npy_intp dims[2];
1153 PyObject *x_i_py = NULL((void*)0);
1154 PyArrayObject *x_i = NULL((void*)0), *BB = NULL((void*)0);
1155 double *t = NULL((void*)0), *h = NULL((void*)0), *ptr;
1156 double x0, xN, arg;
1157
1158 if (!PyArg_ParseTuple(args, "iO", &k, &x_i_py)) {
1159 return NULL((void*)0);
1160 }
1161 if (k < 2) {
1162 PyErr_Format(PyExc_ValueError, "order (%d) must be >=2", k);
1163 return NULL((void*)0);
1164 }
1165
1166 equal = 0;
1167 N = PySequence_LengthPySequence_Size(x_i_py);
1168 if (N == -1 && PyErr_Occurred()) {
1169 PyErr_Clear();
1170 N = PyInt_AsLongPyLong_AsLong(x_i_py);
1171 if (N == -1 && PyErr_Occurred()) {
1172 goto fail;
1173 }
1174 equal = 1;
1175 }
1176 N -= 1;
1177
1178 /* create output matrix */
1179 dims[0] = N + 1;
1180 dims[1] = N + k;
1181 BB = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0)(*(PyObject * (*)(int, npy_intp const *, PyArray_Descr *, int
)) PyArray_API[183])(2, dims, (*(PyArray_Descr * (*)(int)) PyArray_API
[45])(NPY_DOUBLE), 0);
1182 if (BB == NULL((void*)0)) {
1183 goto fail;
1184 }
1185
1186 t = malloc(sizeof(double)*(N + 2*k - 1));
1187 if (t == NULL((void*)0)) {
1188 PyErr_NoMemory();
1189 goto fail;
1190 }
1191
1192 /*
1193 * Create work array to hold computed non-zero values for
1194 * the spline for a value of x.
1195 */
1196 h = malloc(sizeof(double)*(2*k + 1));
1197 if (h == NULL((void*)0)) {
1198 PyErr_NoMemory();
1199 goto fail;
1200 }
1201
1202 numbytes = k*sizeof(double);
1203
1204 if (equal) {
1205 /*
1206 * points equally spaced by 1
1207 * we run deBoor's algorithm one time with artificially created knots
1208 * Then, we keep copying the result to every row
1209 */
1210
1211 /* Create knots at equally-spaced locations from -(K-1) to N+K-1 */
1212 ptr = t;
1213 for (i = -k + 1; i < N + k; i++) {
1214 *ptr++ = i;
1215 }
1216 j = k - 1;
1217 _deBoor_D(t, 0, k, k-1, 0, h);
1218 ptr = PyArray_DATA(BB);
1219 N = N+1;
1220 for (i = 0; i < N; i++) {
1221 memcpy(ptr, h, numbytes);
1222 ptr += N + k;
1223 }
1224 goto finish;
1225 }
1226
1227 /* Not-equally spaced */
1228 x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ARRAY_ALIGNED)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_i_py, (*(PyArray_Descr * (*)
(int)) PyArray_API[45])(NPY_DOUBLE), 1, 1, (((0x0100) & 0x0020
) ? (0x0100) | ((0x0001 | (0x0100 | 0x0400))) : (0x0100)), ((
void*)0));
1229 if (x_i == NULL((void*)0)) {
1230 goto fail;
1231 }
1232 /*
1233 * create dummy knot array with new knots inserted at the end
1234 * selected as mirror symmetric versions of the old knots
1235 */
1236 x0 = *((double *)PyArray_DATA(x_i));
1237 xN = *((double *)PyArray_DATA(x_i) + N);
1238 for (i = 0; i < k - 1; i++) {
1239 /* fill in ends if k > 1*/
1240 t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i, k - 1 - i)((void *)(PyArray_BYTES(x_i) + (k - 1 - i)*PyArray_STRIDES(x_i
)[0]))));
1241 t[k+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i, N - 1 - i)((void *)(PyArray_BYTES(x_i) + (N - 1 - i)*PyArray_STRIDES(x_i
)[0]))));
1242 }
1243 ptr = t + (k - 1);
1244 for (i = 0; i <= N; i++) {
1245 *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)((void *)(PyArray_BYTES(x_i) + (i)*PyArray_STRIDES(x_i)[0]))));
1246 }
1247
1248
1249 /*
1250 * Determine the K+1 non-zero values of the spline and place them in the
1251 * correct location in the matrix for each row (along the diagonals).
1252 * In fact, the last member is always zero so only K non-zero values
1253 * are present.
1254 */
1255 ptr = PyArray_DATA(BB);
1256 for (i = 0, j = k - 1; i < N; i++, j++) {
1257 arg = *((double *)PyArray_DATA(x_i) + i);
1258 _deBoor_D(t, arg, k, j, 0, h);
1259 memcpy(ptr, h, numbytes);
1260 /* advance to next row shifted over one */
1261 ptr += (N + k + 1);
1262 }
1263 /* Last row is different the first coefficient is zero.*/
1264 _deBoor_D(t, xN, k, j - 1, 0, h);
1265 memcpy(ptr, h + 1, numbytes);
1266
1267finish:
1268 Py_XDECREF(x_i)_Py_XDECREF(((PyObject*)(x_i)));
1269 free(t);
1270 free(h);
1271 return (PyObject *)BB;
1272
1273fail:
1274 Py_XDECREF(x_i)_Py_XDECREF(((PyObject*)(x_i)));
1275 Py_XDECREF(BB)_Py_XDECREF(((PyObject*)(BB)));
1276 free(t);
1277 free(h);
1278 return NULL((void*)0);
1279}
1280
1281
1282
1283/*
1284 * Given a set of (N+1) sample positions:
1285 * Construct the (N-1) x (N+K) error matrix J_{ij} such that
1286 *
1287 * for i=1..N-1,
1288 *
1289 * e_i = sum(J_{ij}c_{j},j=-K..N-1)
1290 *
1291 * is the discontinuity of the Kth derivative at the point i in the spline.
1292 *
1293 * This routine takes the N+1 sample positions and the order k and
1294 * constructs the banded matrix J
1295 *
1296 * The returned array is (N+1) times (N+K) ready to be either used
1297 * to compute a minimally Kth-order derivative discontinuous spline
1298 * or to be expanded with an additional K-1 constraints to be used in
1299 * an exact reconstruction approach.
1300 */
1301static char doc_bspldismat[] = "B = _bspldismat(order,xk)\n"
1302"Construct the kth derivative discontinuity jump constraint matrix \n"
1303"for spline fitting of order k given sample positions in xk.\n"
1304"\n"
1305"If xk is an integer (N+1), then the result is equivalent to\n"
1306"xk=arange(N+1)+x0 for any value of x0. This produces the\n"
1307"integer-spaced matrix a bit faster. If xk is a 2-tuple (N+1,dx)\n"
1308"then it produces the result as if the sample distance were dx";
1309static PyObject *
1310_bspldismat(PyObject *dummy, PyObject *args)
1311{
1312 int k, N, i, j, equal, m;
1313 npy_intp dims[2];
1314 PyObject *x_i_py = NULL((void*)0);
1315 PyArrayObject *x_i = NULL((void*)0), *BB = NULL((void*)0);
1316 double *t = NULL((void*)0), *h = NULL((void*)0), *ptr, *dptr;
1317 double x0, xN, dx;
1318
1319 if (!PyArg_ParseTuple(args, "iO", &k, &x_i_py)) {
1320 return NULL((void*)0);
1321 }
1322 if (k < 2) {
1323 PyErr_Format(PyExc_ValueError, "order (%d) must be >=2", k);
1324 return NULL((void*)0);
1325 }
1326
1327 equal = 0;
1328 N = PySequence_LengthPySequence_Size(x_i_py);
1329 if (N == 2 || (N == -1 && PyErr_Occurred())) {
1330 PyErr_Clear();
1331 if (PyTuple_Check(x_i_py)((((((PyObject*)(x_i_py))->ob_type))->tp_flags & ((
1UL << 26))) != 0)) {
1332 /* x_i_py = (N+1, dx) */
1333 N = PyInt_AsLongPyLong_AsLong(PyTuple_GET_ITEM(x_i_py, 0)((((void) (0)), (PyTupleObject *)(x_i_py))->ob_item[0]));
1334 dx = PyFloat_AsDouble(PyTuple_GET_ITEM(x_i_py, 1)((((void) (0)), (PyTupleObject *)(x_i_py))->ob_item[1]));
1335 }
1336 else {
1337 N = PyInt_AsLongPyLong_AsLong(x_i_py);
1338 if (N == -1 && PyErr_Occurred()) {
1339 goto fail;
1340 }
1341 dx = 1.0;
1342 }
1343 equal = 1;
1344 }
1345 N -= 1;
1346
1347 if (N < 2) {
1348 PyErr_Format(PyExc_ValueError, "too few samples (%d)", N);
1349 return NULL((void*)0);
1350 }
1351 /* create output matrix */
1352 dims[0] = N - 1;
1353 dims[1] = N + k;
1354 BB = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0)(*(PyObject * (*)(int, npy_intp const *, PyArray_Descr *, int
)) PyArray_API[183])(2, dims, (*(PyArray_Descr * (*)(int)) PyArray_API
[45])(NPY_DOUBLE), 0);
1355 if (BB == NULL((void*)0)) {
1356 goto fail;
1357 }
1358 t = malloc(sizeof(double)*(N+2*k-1));
1359 if (t == NULL((void*)0)) {
1360 PyErr_NoMemory();
1361 goto fail;
1362 }
1363
1364 /*
1365 * Create work array to hold computed non-zero values for
1366 * the spline for a value of x.
1367 */
1368 h = malloc(sizeof(double)*(2*k + 1));
1369 if (h == NULL((void*)0)) {
1370 PyErr_NoMemory();
1371 goto fail;
1372 }
1373
1374 if (equal) {
1375 /*
1376 * points equally spaced by 1
1377 * we run deBoor's full derivative algorithm twice, subtract the results
1378 * offset by one and then copy the result one time with artificially created knots
1379 * Then, we keep copying the result to every row
1380 */
1381
1382 /* Create knots at equally-spaced locations from -(K-1) to N+K-1 */
1383 double *tmp, factor;
1384 int numbytes;
1385 numbytes = (k + 2)*sizeof(double);
1386 tmp = malloc(numbytes);
1387 if (tmp == NULL((void*)0)) {
1388 PyErr_NoMemory();
1389 goto fail;
1390 }
1391 ptr = t;
1392 for (i = -k + 1; i < N + k; i++) {
1393 *ptr++ = i;
1394 }
1395 j = k - 1;
1396 _deBoor_D(t, 0, k, k-1, k, h);
1397 ptr = tmp;
1398 for (m = 0; m <= k; m++) {
1399 *ptr++ = -h[m];
1400 }
1401 _deBoor_D(t, 0, k, k, k, h);
1402 ptr = tmp + 1;
1403 for (m = 0; m <= k; m++) {
1404 *ptr++ += h[m];
1405 }
1406 if (dx != 1.0) {
1407 factor = pow(dx, (double)k);
1408 for (m = 0; m < k + 2; m++) {
1409 tmp[m] /= factor;
1410 }
1411 }
1412 ptr = PyArray_DATA(BB);
1413 for (i = 0; i < N - 1; i++) {
1414 memcpy(ptr, tmp, numbytes);
1415 ptr += N + k + 1;
1416 }
1417 free(tmp);
1418 goto finish;
1419 }
1420
1421 /* Not-equally spaced */
1422 x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ARRAY_ALIGNED)(*(PyObject * (*)(PyObject *, PyArray_Descr *, int, int, int,
PyObject *)) PyArray_API[69])(x_i_py, (*(PyArray_Descr * (*)
(int)) PyArray_API[45])(NPY_DOUBLE), 1, 1, (((0x0100) & 0x0020
) ? (0x0100) | ((0x0001 | (0x0100 | 0x0400))) : (0x0100)), ((
void*)0));
1423 if (x_i == NULL((void*)0)) {
1424 goto fail;
1425 }
1426 /*
1427 * create dummy knot array with new knots inserted at the end
1428 * selected as mirror symmetric versions of the old knots
1429 */
1430 x0 = *((double *)PyArray_DATA(x_i));
1431 xN = *((double *)PyArray_DATA(x_i) + N);
1432 for (i = 0; i < k - 1; i++) {
1433 /* fill in ends if k > 1*/
1434 t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i, k - 1 - i)((void *)(PyArray_BYTES(x_i) + (k - 1 - i)*PyArray_STRIDES(x_i
)[0]))));
1435 t[k+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i, N - 1 - i)((void *)(PyArray_BYTES(x_i) + (N - 1 - i)*PyArray_STRIDES(x_i
)[0]))));
1436 }
1437 ptr = t + (k - 1);
1438 for (i = 0; i <= N; i++) {
1439 *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)((void *)(PyArray_BYTES(x_i) + (i)*PyArray_STRIDES(x_i)[0]))));
1440 }
1441
1442
1443 /*
1444 * Determine the K+1 non-zero values of the discontinuity jump matrix
1445 * and place them in the correct location in the matrix for each row
1446 * (along the diagonals).
1447 *
1448 * The matrix is
1449 *
1450 * J_{ij} = b^{(k)}_{j,k}(x^{+}_i) - b^{(k)}_{j,k}(x^{-}_i)
1451 */
1452 ptr = PyArray_DATA(BB);
1453 dptr = ptr;
1454 for (i = 0, j = k - 1; i < N - 1; i++, j++) {
1455 _deBoor_D(t, 0, k, j, k, h);
1456 /* We need to copy over but negate the terms */
1457 for (m = 0; m <= k; m++) {
1458 *ptr++ = -h[m];
1459 }
1460 /*
1461 * If we are past the first row, then we need to also add the current
1462 * values result to the previous row
1463 */
1464 if (i > 0) {
1465 for (m = 0; m <= k; m++) {
1466 *dptr++ += h[m];
1467 }
1468 }
1469 /* store location of last start position plus one.*/
1470 dptr = ptr - k;
1471 /* advance to next row shifted over one */
1472 ptr += N;
1473 }
1474 /* We need to finish the result for the last row. */
1475 _deBoor_D(t, 0, k, j, k, h);
1476 for (m = 0; m <= k; m++) {
1477 *dptr++ += h[m];
1478 }
1479
1480finish:
1481 Py_XDECREF(x_i)_Py_XDECREF(((PyObject*)(x_i)));
1482 free(t);
1483 free(h);
1484 return (PyObject *)BB;
1485
1486fail:
1487 Py_XDECREF(x_i)_Py_XDECREF(((PyObject*)(x_i)));
1488 Py_XDECREF(BB)_Py_XDECREF(((PyObject*)(BB)));
1489 free(t);
1490 free(h);
1491 return NULL((void*)0);
1492}
1493
1494/* End of functions moved verbatim from __fitpack.h */
1495
1496
1497
1498static struct PyMethodDef fitpack_module_methods[] = {
1499{"_curfit",
1500 fitpack_curfit,
1501 METH_VARARGS0x0001, doc_curfit},
1502{"_spl_",
1503 fitpack_spl_,
1504 METH_VARARGS0x0001, doc_spl_},
1505{"_splint",
1506 fitpack_splint,
1507 METH_VARARGS0x0001, doc_splint},
1508{"_sproot",
1509 fitpack_sproot,
1510 METH_VARARGS0x0001, doc_sproot},
1511{"_spalde",
1512 fitpack_spalde,
1513 METH_VARARGS0x0001, doc_spalde},
1514{"_parcur",
1515 fitpack_parcur,
1516 METH_VARARGS0x0001, doc_parcur},
1517{"_surfit",
1518 fitpack_surfit,
1519 METH_VARARGS0x0001, doc_surfit},
1520{"_bispev",
1521 fitpack_bispev,
1522 METH_VARARGS0x0001, doc_bispev},
1523{"_insert",
1524 fitpack_insert,
1525 METH_VARARGS0x0001, doc_insert},
1526{"_bspleval",
1527 _bspleval,
1528 METH_VARARGS0x0001, doc_bspleval},
1529{"_bsplmat",
1530 _bsplmat,
1531 METH_VARARGS0x0001, doc_bsplmat},
1532{"_bspldismat",
1533 _bspldismat,
1534 METH_VARARGS0x0001, doc_bspldismat},
1535{NULL((void*)0), NULL((void*)0), 0, NULL((void*)0)}
1536};
1537
1538static struct PyModuleDef moduledef = {
1539 PyModuleDef_HEAD_INIT{ { 1, ((void*)0) }, ((void*)0), 0, ((void*)0), },
1540 "_fitpack",
1541 NULL((void*)0),
1542 -1,
1543 fitpack_module_methods,
1544 NULL((void*)0),
1545 NULL((void*)0),
1546 NULL((void*)0),
1547 NULL((void*)0)
1548};
1549
1550PyObject *PyInit__fitpack(void)
1551{
1552 PyObject *m, *d, *s;
1553
1554 m = PyModule_Create(&moduledef)PyModule_Create2(&moduledef, 1013);
1
Calling 'PyModule_Create2'
→
3
←
Returning from 'PyModule_Create2'
→
5
←
PyObject ownership leak with reference count of 1
1555 import_array(){if (_import_array() < 0) {PyErr_Print(); PyErr_SetString(
PyExc_ImportError, "numpy.core.multiarray failed to import");
return ((void*)0); } };
4
←
Taking true branch
→
1556
1557 d = PyModule_GetDict(m);
1558
1559 s = PyUnicode_FromString(" 1.7 ");
1560 PyDict_SetItemString(d, "__version__", s);
1561 fitpack_error = PyErr_NewException ("fitpack.error", NULL((void*)0), NULL((void*)0));
1562 Py_DECREF(s)_Py_DECREF(((PyObject*)(s)));
1563 if (PyErr_Occurred()) {
1564 Py_FatalError("can't initialize module fitpack");
1565 }
1566
1567 return m;
1568}
←

/opt/pyrefcon/lib/pyrefcon/models/models/PyModule_Create2.model

1#ifndef PyModule_Create2
2struct _object;
3typedef struct _object PyObject;
4PyObject* clang_analyzer_PyObject_New_Reference();
5PyObject* PyModule_Create2(PyModuleDef *def, int module_api_version) {
6 return clang_analyzer_PyObject_New_Reference();
2
←
Setting reference count to 1
→
7}
8#else
9#warning "API PyModule_Create2 is defined as a macro."
10#endif