1 /* This file is used to make _multipackmodule.c */
2 /* $Revision$ */
3 /* module_methods:
4 {"_hybrd", minpack_hybrd, METH_VARARGS, doc_hybrd},
5 {"_hybrj", minpack_hybrj, METH_VARARGS, doc_hybrj},
6 {"_lmdif", minpack_lmdif, METH_VARARGS, doc_lmdif},
7 {"_lmder", minpack_lmder, METH_VARARGS, doc_lmder},
8 {"_chkder", minpack_chkder, METH_VARARGS, doc_chkder},
9 */
10
11 /* link libraries:
12 minpack
13 linpack_lite
14 blas
15 */
16
17 /* python files:
18 minpack.py
19 */
20
21 #if defined(NO_APPEND_FORTRAN)
22 #if defined(UPPERCASE_FORTRAN)
23 /* nothing to do in that case */
24 #else
25 #define CHKDER chkder
26 #define HYBRD hybrd
27 #define HYBRJ hybrj
28 #define LMDIF lmdif
29 #define LMDER lmder
30 #define LMSTR lmstr
31 #endif
32 #else
33 #if defined(UPPERCASE_FORTRAN)
34 #define CHKDER CHKDER_
35 #define HYBRD HYBRD_
36 #define HYBRJ HYBRJ_
37 #define LMDIF LMDIF_
38 #define LMDER LMDER_
39 #define LMSTR LMSTR_
40 #else
41 #define CHKDER chkder_
42 #define HYBRD hybrd_
43 #define HYBRJ hybrj_
44 #define LMDIF lmdif_
45 #define LMDER lmder_
46 #define LMSTR lmstr_
47 #endif
48 #endif
49
50 extern void CHKDER(int*,int*,double*,double*,double*,int*,double*,double*,int*,double*);
51 extern void HYBRD(void*,int*,double*,double*,double*,int*,int*,int*,double*,double*,int*,double*,int*,int*,int*,double*,int*,double*,int*,double*,double*,double*,double*,double*);
52 extern void HYBRJ(void*,int*,double*,double*,double*,int*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,int*,double*,double*,double*,double*,double*);
53 extern void LMDIF(void*,int*,int*,double*,double*,double*,double*,double*,int*,double*,double*,int*,double*,int*,int*,int*,double*,int*,int*,double*,double*,double*,double*,double*);
54 extern void LMDER(void*,int*,int*,double*,double*,double*,int*,double*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,int*,double*,double*,double*,double*,double*);
55 extern void LMSTR(void*,int*,int*,double*,double*,double*,int*,double*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,int*,double*,double*,double*,double*,double*);
56
57 /* We only use ccallback with Python functions right now */
58 static ccallback_signature_t call_signatures[] = {
59 {NULL}
60 };
61
init_callback(ccallback_t * callback,PyObject * fcn,PyObject * extra_args)62 static int init_callback(ccallback_t *callback, PyObject *fcn, PyObject *extra_args)
63 {
64 int ret;
65 int flags = CCALLBACK_OBTAIN;
66
67 ret = ccallback_prepare(callback, call_signatures, fcn, flags);
68 if (ret == -1) {
69 return -1;
70 }
71
72 callback->info_p = (void *)extra_args;
73
74 return 0;
75 }
76
release_callback(ccallback_t * callback)77 static int release_callback(ccallback_t *callback)
78 {
79 return ccallback_release(callback) != 0;
80 }
81
init_jac_callback(ccallback_t * callback,jac_callback_info_t * jac_callback_info,PyObject * fcn,PyObject * Dfun,PyObject * extra_args,int col_deriv)82 static int init_jac_callback(ccallback_t *callback, jac_callback_info_t *jac_callback_info, PyObject *fcn, PyObject *Dfun, PyObject *extra_args, int col_deriv)
83 {
84 int ret;
85 int flags = CCALLBACK_OBTAIN;
86
87 ret = ccallback_prepare(callback, call_signatures, fcn, flags);
88 if (ret == -1) {
89 return -1;
90 }
91
92 jac_callback_info->Dfun = Dfun;
93 jac_callback_info->extra_args = extra_args;
94 jac_callback_info->jac_transpose = !col_deriv;
95
96 callback->info_p = (void *)jac_callback_info;
97
98 return 0;
99 }
100
101
raw_multipack_calling_function(int * n,double * x,double * fvec,int * iflag)102 int raw_multipack_calling_function(int *n, double *x, double *fvec, int *iflag)
103 {
104 /* This is the function called from the Fortran code it should
105 -- use call_python_function to get a multiarrayobject result
106 -- check for errors and return -1 if any
107 -- otherwise place result of calculation in *fvec
108 */
109
110 ccallback_t *callback = ccallback_obtain();
111 PyObject *multipack_python_function = callback->py_function;
112 PyObject *multipack_extra_arguments = (PyObject *)callback->info_p;
113
114 PyArrayObject *result_array = NULL;
115
116 result_array = (PyArrayObject *)call_python_function(multipack_python_function, *n, x, multipack_extra_arguments, 1, minpack_error, *n);
117 if (result_array == NULL) {
118 *iflag = -1;
119 return -1;
120 }
121 memcpy(fvec, PyArray_DATA(result_array), (*n)*sizeof(double));
122 Py_DECREF(result_array);
123 return 0;
124
125 }
126
127
jac_multipack_calling_function(int * n,double * x,double * fvec,double * fjac,int * ldfjac,int * iflag)128 int jac_multipack_calling_function(int *n, double *x, double *fvec, double *fjac, int *ldfjac, int *iflag)
129 {
130 /* This is the function called from the Fortran code it should
131 -- use call_python_function to get a multiarrayobject result
132 -- check for errors and return -1 if any
133 -- otherwise place result of calculation in *fvec or *fjac.
134
135 If iflag = 1 this should compute the function.
136 If iflag = 2 this should compute the jacobian (derivative matrix)
137 */
138
139 ccallback_t *callback = ccallback_obtain();
140 PyObject *multipack_python_function = callback->py_function,
141 *multipack_python_jacobian = ((jac_callback_info_t *)callback->info_p)->Dfun;
142 PyObject *multipack_extra_arguments = ((jac_callback_info_t *)callback->info_p)->extra_args;
143 int multipack_jac_transpose = ((jac_callback_info_t *)callback->info_p)->jac_transpose;
144
145 PyArrayObject *result_array;
146
147 if (*iflag == 1) {
148 result_array = (PyArrayObject *)call_python_function(multipack_python_function, *n, x, multipack_extra_arguments, 1, minpack_error, *n);
149 if (result_array == NULL) {
150 *iflag = -1;
151 return -1;
152 }
153 memcpy(fvec, PyArray_DATA(result_array), (*n)*sizeof(double));
154 }
155 else { /* iflag == 2 */
156 result_array = (PyArrayObject *)call_python_function(multipack_python_jacobian, *n, x, multipack_extra_arguments, 2, minpack_error, (*n)*(*ldfjac));
157 if (result_array == NULL) {
158 *iflag = -1;
159 return -1;
160 }
161 if (multipack_jac_transpose == 1)
162 MATRIXC2F(fjac, PyArray_DATA(result_array), *n, *ldfjac)
163 else
164 memcpy(fjac, PyArray_DATA(result_array), (*n)*(*ldfjac)*sizeof(double));
165 }
166
167 Py_DECREF(result_array);
168 return 0;
169 }
170
raw_multipack_lm_function(int * m,int * n,double * x,double * fvec,int * iflag)171 int raw_multipack_lm_function(int *m, int *n, double *x, double *fvec, int *iflag)
172 {
173 /* This is the function called from the Fortran code it should
174 -- use call_python_function to get a multiarrayobject result
175 -- check for errors and return -1 if any
176 -- otherwise place result of calculation in *fvec
177 */
178
179 ccallback_t *callback = ccallback_obtain();
180 PyObject *multipack_python_function = callback->py_function;
181 PyObject *multipack_extra_arguments = (PyObject *)callback->info_p;
182
183 PyArrayObject *result_array = NULL;
184
185 result_array = (PyArrayObject *)call_python_function(multipack_python_function,*n, x, multipack_extra_arguments, 1, minpack_error, *m);
186 if (result_array == NULL) {
187 *iflag = -1;
188 return -1;
189 }
190 memcpy(fvec, PyArray_DATA(result_array), (*m)*sizeof(double));
191 Py_DECREF(result_array);
192 return 0;
193 }
194
jac_multipack_lm_function(int * m,int * n,double * x,double * fvec,double * fjac,int * ldfjac,int * iflag)195 int jac_multipack_lm_function(int *m, int *n, double *x, double *fvec, double *fjac, int *ldfjac, int *iflag)
196 {
197 /* This is the function called from the Fortran code it should
198 -- use call_python_function to get a multiarrayobject result
199 -- check for errors and return -1 if any
200 -- otherwise place result of calculation in *fvec or *fjac.
201
202 If iflag = 1 this should compute the function.
203 If iflag = 2 this should compute the jacobian (derivative matrix)
204 */
205
206 ccallback_t *callback = ccallback_obtain();
207 PyObject *multipack_python_function = callback->py_function,
208 *multipack_python_jacobian = ((jac_callback_info_t *)callback->info_p)->Dfun;
209 PyObject *multipack_extra_arguments = ((jac_callback_info_t *)callback->info_p)->extra_args;
210 int multipack_jac_transpose = ((jac_callback_info_t *)callback->info_p)->jac_transpose;
211
212 PyArrayObject *result_array;
213
214 if (*iflag == 1) {
215 result_array = (PyArrayObject *)call_python_function(multipack_python_function, *n, x, multipack_extra_arguments, 1, minpack_error, *m);
216 if (result_array == NULL) {
217 *iflag = -1;
218 return -1;
219 }
220 memcpy(fvec, PyArray_DATA(result_array), (*m)*sizeof(double));
221 }
222 else { /* iflag == 2 */
223 result_array = (PyArrayObject *)call_python_function(multipack_python_jacobian, *n, x, multipack_extra_arguments, 2, minpack_error, (*n)*(*ldfjac));
224 if (result_array == NULL) {
225 *iflag = -1;
226 return -1;
227 }
228 if (multipack_jac_transpose == 1)
229 MATRIXC2F(fjac, PyArray_DATA(result_array), *n, *ldfjac)
230 else
231 memcpy(fjac, PyArray_DATA(result_array), (*n)*(*ldfjac)*sizeof(double));
232 }
233
234 Py_DECREF(result_array);
235 return 0;
236 }
237
238
239 static char doc_hybrd[] = "[x,infodict,info] = _hybrd(fun, x0, args, full_output, xtol, maxfev, ml, mu, epsfcn, factor, diag)";
240
minpack_hybrd(PyObject * dummy,PyObject * args)241 static PyObject *minpack_hybrd(PyObject *dummy, PyObject *args) {
242 PyObject *fcn, *x0, *extra_args = NULL, *o_diag = NULL;
243 int full_output = 0, maxfev = -10, ml = -10, mu = -10;
244 double xtol = 1.49012e-8, epsfcn = 0.0, factor = 1.0e2;
245 int mode = 2, nprint = 0, info, nfev, ldfjac;
246 npy_intp n,lr;
247 int n_int, lr_int; /* for casted storage to pass int into HYBRD */
248 double *x, *fvec, *diag, *fjac, *r, *qtf;
249
250 PyArrayObject *ap_x = NULL, *ap_fvec = NULL;
251 PyArrayObject *ap_fjac = NULL, *ap_r = NULL, *ap_qtf = NULL;
252 PyArrayObject *ap_diag = NULL;
253
254 npy_intp dims[2];
255 int allocated = 0;
256 double *wa = NULL;
257
258 STORE_VARS(); /* Define storage variables for global variables. */
259
260 if (!PyArg_ParseTuple(args, "OO|OidiiiddO", &fcn, &x0, &extra_args, &full_output, &xtol, &maxfev, &ml, &mu, &epsfcn, &factor, &o_diag)) return NULL;
261
262 INIT_FUNC(fcn,extra_args,minpack_error);
263
264 /* Initial input vector */
265 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x0, NPY_DOUBLE, 1, 1);
266 if (ap_x == NULL) goto fail;
267 x = (double *) PyArray_DATA(ap_x);
268 n = PyArray_DIMS(ap_x)[0];
269
270 lr = n * (n + 1) / 2;
271 if (ml < 0) ml = n-1;
272 if (mu < 0) mu = n-1;
273 if (maxfev < 0) maxfev = 200*(n+1);
274
275 /* Setup array to hold the function evaluations */
276 ap_fvec = (PyArrayObject *)call_python_function(fcn, n, x, extra_args, 1, minpack_error, -1);
277 if (ap_fvec == NULL) goto fail;
278 fvec = (double *) PyArray_DATA(ap_fvec);
279 if (PyArray_NDIM(ap_fvec) == 0)
280 n = 1;
281 else if (PyArray_DIMS(ap_fvec)[0] < n)
282 n = PyArray_DIMS(ap_fvec)[0];
283
284 SET_DIAG(ap_diag,o_diag,mode);
285
286 dims[0] = n; dims[1] = n;
287 ap_r = (PyArrayObject *)PyArray_SimpleNew(1,&lr,NPY_DOUBLE);
288 ap_qtf = (PyArrayObject *)PyArray_SimpleNew(1,&n,NPY_DOUBLE);
289 ap_fjac = (PyArrayObject *)PyArray_SimpleNew(2,dims,NPY_DOUBLE);
290
291 if (ap_r == NULL || ap_qtf == NULL || ap_fjac ==NULL) goto fail;
292
293 r = (double *) PyArray_DATA(ap_r);
294 qtf = (double *) PyArray_DATA(ap_qtf);
295 fjac = (double *) PyArray_DATA(ap_fjac);
296 ldfjac = dims[1];
297
298 if ((wa = malloc(4*n * sizeof(double)))==NULL) {
299 PyErr_NoMemory();
300 goto fail;
301 }
302 allocated = 1;
303
304 /* Call the underlying FORTRAN routines. */
305 n_int = n; lr_int = lr; /* cast/store/pass into HYBRD */
306 HYBRD(raw_multipack_calling_function, &n_int, x, fvec, &xtol, &maxfev, &ml, &mu, &epsfcn, diag, &mode, &factor, &nprint, &info, &nfev, fjac, &ldfjac, r, &lr_int, qtf, wa, wa+n, wa+2*n, wa+3*n);
307
308 RESTORE_FUNC();
309
310 if (info < 0) goto fail; /* Python Terminated */
311
312
313 free(wa);
314 Py_DECREF(extra_args);
315 Py_DECREF(ap_diag);
316
317 if (full_output) {
318 return Py_BuildValue("N{s:N,s:i,s:N,s:N,s:N}i",PyArray_Return(ap_x),"fvec",PyArray_Return(ap_fvec),"nfev",nfev,"fjac",PyArray_Return(ap_fjac),"r",PyArray_Return(ap_r),"qtf",PyArray_Return(ap_qtf),info);
319 }
320 else {
321 Py_DECREF(ap_fvec);
322 Py_DECREF(ap_fjac);
323 Py_DECREF(ap_r);
324 Py_DECREF(ap_qtf);
325 return Py_BuildValue("Ni",PyArray_Return(ap_x),info);
326 }
327
328 fail:
329 RESTORE_FUNC();
330 fail_free:
331 Py_XDECREF(extra_args);
332 Py_XDECREF(ap_x);
333 Py_XDECREF(ap_fvec);
334 Py_XDECREF(ap_diag);
335 Py_XDECREF(ap_fjac);
336 Py_XDECREF(ap_r);
337 Py_XDECREF(ap_qtf);
338 if (allocated) free(wa);
339 return NULL;
340 }
341
342
343 static char doc_hybrj[] = "[x,infodict,info] = _hybrj(fun, Dfun, x0, args, full_output, col_deriv, xtol, maxfev, factor, diag)";
344
minpack_hybrj(PyObject * dummy,PyObject * args)345 static PyObject *minpack_hybrj(PyObject *dummy, PyObject *args) {
346 PyObject *fcn, *Dfun, *x0, *extra_args = NULL, *o_diag = NULL;
347 int full_output = 0, maxfev = -10, col_deriv = 1;
348 double xtol = 1.49012e-8, factor = 1.0e2;
349 int mode = 2, nprint = 0, info, nfev, njev, ldfjac;
350 npy_intp n, lr;
351 int n_int, lr_int;
352 double *x, *fvec, *diag, *fjac, *r, *qtf;
353
354 PyArrayObject *ap_x = NULL, *ap_fvec = NULL;
355 PyArrayObject *ap_fjac = NULL, *ap_r = NULL, *ap_qtf = NULL;
356 PyArrayObject *ap_diag = NULL;
357
358 npy_intp dims[2];
359 int allocated = 0;
360 double *wa = NULL;
361
362 STORE_VARS();
363
364 if (!PyArg_ParseTuple(args, "OOO|OiididO", &fcn, &Dfun, &x0, &extra_args, &full_output, &col_deriv, &xtol, &maxfev, &factor, &o_diag)) return NULL;
365
366 INIT_JAC_FUNC(fcn,Dfun,extra_args,col_deriv,minpack_error);
367
368 /* Initial input vector */
369 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x0, NPY_DOUBLE, 1, 1);
370 if (ap_x == NULL) goto fail;
371 x = (double *) PyArray_DATA(ap_x);
372 n = PyArray_DIMS(ap_x)[0];
373 lr = n * (n + 1) / 2;
374
375 if (maxfev < 0) maxfev = 100*(n+1);
376
377 /* Setup array to hold the function evaluations */
378 ap_fvec = (PyArrayObject *)call_python_function(fcn, n, x, extra_args, 1, minpack_error, -1);
379 if (ap_fvec == NULL) goto fail;
380 fvec = (double *) PyArray_DATA(ap_fvec);
381 if (PyArray_NDIM(ap_fvec) == 0)
382 n = 1;
383 else if (PyArray_DIMS(ap_fvec)[0] < n)
384 n = PyArray_DIMS(ap_fvec)[0];
385
386 SET_DIAG(ap_diag,o_diag,mode);
387
388 dims[0] = n; dims[1] = n;
389 ap_r = (PyArrayObject *)PyArray_SimpleNew(1,&lr,NPY_DOUBLE);
390 ap_qtf = (PyArrayObject *)PyArray_SimpleNew(1,&n,NPY_DOUBLE);
391 ap_fjac = (PyArrayObject *)PyArray_SimpleNew(2,dims,NPY_DOUBLE);
392
393 if (ap_r == NULL || ap_qtf == NULL || ap_fjac ==NULL) goto fail;
394
395 r = (double *) PyArray_DATA(ap_r);
396 qtf = (double *) PyArray_DATA(ap_qtf);
397 fjac = (double *) PyArray_DATA(ap_fjac);
398
399 ldfjac = dims[1];
400
401 if ((wa = malloc(4*n * sizeof(double)))==NULL) {
402 PyErr_NoMemory();
403 goto fail;
404 }
405 allocated = 1;
406
407 /* Call the underlying FORTRAN routines. */
408 n_int = n; lr_int = lr; /* cast/store/pass into HYBRJ */
409 HYBRJ(jac_multipack_calling_function, &n_int, x, fvec, fjac, &ldfjac, &xtol, &maxfev, diag, &mode, &factor, &nprint, &info, &nfev, &njev, r, &lr_int, qtf, wa, wa+n, wa+2*n, wa+3*n);
410
411 RESTORE_JAC_FUNC();
412
413 if (info < 0) goto fail; /* Python Terminated */
414
415 free(wa);
416 Py_DECREF(extra_args);
417 Py_DECREF(ap_diag);
418
419 if (full_output) {
420 return Py_BuildValue("N{s:N,s:i,s:i,s:N,s:N,s:N}i",PyArray_Return(ap_x),"fvec",PyArray_Return(ap_fvec),"nfev",nfev,"njev",njev,"fjac",PyArray_Return(ap_fjac),"r",PyArray_Return(ap_r),"qtf",PyArray_Return(ap_qtf),info);
421 }
422 else {
423 Py_DECREF(ap_fvec);
424 Py_DECREF(ap_fjac);
425 Py_DECREF(ap_r);
426 Py_DECREF(ap_qtf);
427 return Py_BuildValue("Ni",PyArray_Return(ap_x),info);
428 }
429
430 fail:
431 RESTORE_JAC_FUNC();
432 fail_free:
433 Py_XDECREF(extra_args);
434 Py_XDECREF(ap_x);
435 Py_XDECREF(ap_fvec);
436 Py_XDECREF(ap_fjac);
437 Py_XDECREF(ap_diag);
438 Py_XDECREF(ap_r);
439 Py_XDECREF(ap_qtf);
440 if (allocated) free(wa);
441 return NULL;
442
443 }
444
445 /************************ Levenberg-Marquardt *******************/
446
447 static char doc_lmdif[] = "[x,infodict,info] = _lmdif(fun, x0, args, full_output, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)";
448
minpack_lmdif(PyObject * dummy,PyObject * args)449 static PyObject *minpack_lmdif(PyObject *dummy, PyObject *args) {
450 PyObject *fcn, *x0, *extra_args = NULL, *o_diag = NULL;
451 int full_output = 0, maxfev = -10;
452 double xtol = 1.49012e-8, ftol = 1.49012e-8;
453 double gtol = 0.0, epsfcn = 0.0, factor = 1.0e2;
454 int m, mode = 2, nprint = 0, info = 0, nfev, ldfjac, *ipvt;
455 npy_intp n;
456 int n_int; /* for casted storage to pass int into LMDIF */
457 double *x, *fvec, *diag, *fjac, *qtf;
458
459 PyArrayObject *ap_x = NULL, *ap_fvec = NULL;
460 PyArrayObject *ap_fjac = NULL, *ap_ipvt = NULL, *ap_qtf = NULL;
461 PyArrayObject *ap_diag = NULL;
462
463 npy_intp dims[2];
464 int allocated = 0;
465 double *wa = NULL;
466
467 STORE_VARS();
468
469 if (!PyArg_ParseTuple(args, "OO|OidddiddO", &fcn, &x0, &extra_args, &full_output, &ftol, &xtol, >ol, &maxfev, &epsfcn, &factor, &o_diag)) return NULL;
470
471 INIT_FUNC(fcn,extra_args,minpack_error);
472
473 /* Initial input vector */
474 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x0, NPY_DOUBLE, 1, 1);
475 if (ap_x == NULL) goto fail;
476 x = (double *) PyArray_DATA(ap_x);
477 n = PyArray_DIMS(ap_x)[0];
478 dims[0] = n;
479
480 SET_DIAG(ap_diag,o_diag,mode);
481
482 if (maxfev < 0) maxfev = 200*(n+1);
483
484 /* Setup array to hold the function evaluations and find it's size*/
485 ap_fvec = (PyArrayObject *)call_python_function(fcn, n, x, extra_args, 1, minpack_error, -1);
486 if (ap_fvec == NULL) goto fail;
487 fvec = (double *) PyArray_DATA(ap_fvec);
488 m = (PyArray_NDIM(ap_fvec) > 0 ? PyArray_DIMS(ap_fvec)[0] : 1);
489
490 dims[0] = n; dims[1] = m;
491 ap_ipvt = (PyArrayObject *)PyArray_SimpleNew(1,&n,NPY_INT);
492 ap_qtf = (PyArrayObject *)PyArray_SimpleNew(1,&n,NPY_DOUBLE);
493 ap_fjac = (PyArrayObject *)PyArray_SimpleNew(2,dims,NPY_DOUBLE);
494
495 if (ap_ipvt == NULL || ap_qtf == NULL || ap_fjac ==NULL) goto fail;
496
497 ipvt = (int *) PyArray_DATA(ap_ipvt);
498 qtf = (double *) PyArray_DATA(ap_qtf);
499 fjac = (double *) PyArray_DATA(ap_fjac);
500 ldfjac = dims[1];
501 wa = (double *)malloc((3*n + m)* sizeof(double));
502 if (wa == NULL) {
503 PyErr_NoMemory();
504 goto fail;
505 }
506 allocated = 1;
507
508 /* Call the underlying FORTRAN routines. */
509 n_int = n; /* to provide int*-pointed storage for int argument of LMDIF */
510 LMDIF(raw_multipack_lm_function, &m, &n_int, x, fvec, &ftol, &xtol, >ol, &maxfev, &epsfcn, diag, &mode, &factor, &nprint, &info, &nfev, fjac, &ldfjac, ipvt, qtf, wa, wa+n, wa+2*n, wa+3*n);
511
512 RESTORE_FUNC();
513
514 if (info < 0) goto fail; /* Python error */
515
516 free(wa);
517 Py_DECREF(extra_args);
518 Py_DECREF(ap_diag);
519
520 if (full_output) {
521 return Py_BuildValue("N{s:N,s:i,s:N,s:N,s:N}i",PyArray_Return(ap_x),"fvec",PyArray_Return(ap_fvec),"nfev",nfev,"fjac",PyArray_Return(ap_fjac),"ipvt",PyArray_Return(ap_ipvt),"qtf",PyArray_Return(ap_qtf),info);
522 }
523 else {
524 Py_DECREF(ap_fvec);
525 Py_DECREF(ap_fjac);
526 Py_DECREF(ap_ipvt);
527 Py_DECREF(ap_qtf);
528 return Py_BuildValue("Ni",PyArray_Return(ap_x),info);
529 }
530
531 fail:
532 RESTORE_FUNC();
533 fail_free:
534 Py_XDECREF(extra_args);
535 Py_XDECREF(ap_x);
536 Py_XDECREF(ap_fvec);
537 Py_XDECREF(ap_fjac);
538 Py_XDECREF(ap_diag);
539 Py_XDECREF(ap_ipvt);
540 Py_XDECREF(ap_qtf);
541 if (allocated) free(wa);
542 return NULL;
543 }
544
545
546 static char doc_lmder[] = "[x,infodict,info] = _lmder(fun, Dfun, x0, args, full_output, col_deriv, ftol, xtol, gtol, maxfev, factor, diag)";
547
minpack_lmder(PyObject * dummy,PyObject * args)548 static PyObject *minpack_lmder(PyObject *dummy, PyObject *args) {
549 PyObject *fcn, *x0, *Dfun, *extra_args = NULL, *o_diag = NULL;
550 int full_output = 0, maxfev = -10, col_deriv = 1;
551 double xtol = 1.49012e-8, ftol = 1.49012e-8;
552 double gtol = 0.0, factor = 1.0e2;
553 int m, mode = 2, nprint = 0, info, nfev, njev, ldfjac, *ipvt;
554 npy_intp n;
555 int n_int;
556 double *x, *fvec, *diag, *fjac, *qtf;
557
558 PyArrayObject *ap_x = NULL, *ap_fvec = NULL;
559 PyArrayObject *ap_fjac = NULL, *ap_ipvt = NULL, *ap_qtf = NULL;
560 PyArrayObject *ap_diag = NULL;
561
562 npy_intp dims[2];
563 int allocated = 0;
564 double *wa = NULL;
565
566 STORE_VARS();
567
568 if (!PyArg_ParseTuple(args, "OOO|OiidddidO", &fcn, &Dfun, &x0, &extra_args, &full_output, &col_deriv, &ftol, &xtol, >ol, &maxfev, &factor, &o_diag)) return NULL;
569
570 INIT_JAC_FUNC(fcn,Dfun,extra_args,col_deriv,minpack_error);
571
572 /* Initial input vector */
573 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x0, NPY_DOUBLE, 1, 1);
574 if (ap_x == NULL) goto fail;
575 x = (double *) PyArray_DATA(ap_x);
576 n = PyArray_DIMS(ap_x)[0];
577
578 if (maxfev < 0) maxfev = 100*(n+1);
579
580 /* Setup array to hold the function evaluations */
581 ap_fvec = (PyArrayObject *)call_python_function(fcn, n, x, extra_args, 1, minpack_error, -1);
582 if (ap_fvec == NULL) goto fail;
583 fvec = (double *) PyArray_DATA(ap_fvec);
584
585 SET_DIAG(ap_diag,o_diag,mode);
586
587 m = (PyArray_NDIM(ap_fvec) > 0 ? PyArray_DIMS(ap_fvec)[0] : 1);
588
589 dims[0] = n; dims[1] = m;
590 ap_ipvt = (PyArrayObject *)PyArray_SimpleNew(1,&n,NPY_INT);
591 ap_qtf = (PyArrayObject *)PyArray_SimpleNew(1,&n,NPY_DOUBLE);
592 ap_fjac = (PyArrayObject *)PyArray_SimpleNew(2,dims,NPY_DOUBLE);
593
594 if (ap_ipvt == NULL || ap_qtf == NULL || ap_fjac ==NULL) goto fail;
595
596 ipvt = (int *) PyArray_DATA(ap_ipvt);
597 qtf = (double *) PyArray_DATA(ap_qtf);
598 fjac = (double *) PyArray_DATA(ap_fjac);
599 ldfjac = dims[1];
600 wa = (double *)malloc((3*n + m)* sizeof(double));
601 if (wa == NULL) {
602 PyErr_NoMemory();
603 goto fail;
604 }
605 allocated = 1;
606
607 /* Call the underlying FORTRAN routines. */
608 n_int = n;
609 LMDER(jac_multipack_lm_function, &m, &n_int, x, fvec, fjac, &ldfjac, &ftol, &xtol, >ol, &maxfev, diag, &mode, &factor, &nprint, &info, &nfev, &njev, ipvt, qtf, wa, wa+n, wa+2*n, wa+3*n);
610
611 RESTORE_JAC_FUNC();
612
613 if (info < 0) goto fail; /* Python error */
614
615 free(wa);
616 Py_DECREF(extra_args);
617 Py_DECREF(ap_diag);
618
619 if (full_output) {
620 return Py_BuildValue("N{s:N,s:i,s:i,s:N,s:N,s:N}i",PyArray_Return(ap_x),"fvec",PyArray_Return(ap_fvec),"nfev",nfev,"njev",njev,"fjac",PyArray_Return(ap_fjac),"ipvt",PyArray_Return(ap_ipvt),"qtf",PyArray_Return(ap_qtf),info);
621 }
622 else {
623 Py_DECREF(ap_fvec);
624 Py_DECREF(ap_fjac);
625 Py_DECREF(ap_ipvt);
626 Py_DECREF(ap_qtf);
627 return Py_BuildValue("Ni",PyArray_Return(ap_x),info);
628 }
629
630 fail:
631 RESTORE_JAC_FUNC();
632 fail_free:
633 Py_XDECREF(extra_args);
634 Py_XDECREF(ap_x);
635 Py_XDECREF(ap_fvec);
636 Py_XDECREF(ap_fjac);
637 Py_XDECREF(ap_diag);
638 Py_XDECREF(ap_ipvt);
639 Py_XDECREF(ap_qtf);
640 if (allocated) free(wa);
641 return NULL;
642 }
643
644
645 /** Check gradient function **/
646
647 static char doc_chkder[] = "_chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err)";
648
minpack_chkder(PyObject * self,PyObject * args)649 static PyObject *minpack_chkder(PyObject *self, PyObject *args)
650 {
651 PyArrayObject *ap_fvecp = NULL, *ap_fjac = NULL, *ap_err = NULL;
652 PyArrayObject *ap_x = NULL, *ap_fvec = NULL, *ap_xp = NULL;
653 PyObject *o_x, *o_fvec, *o_fjac, *o_fvecp;
654 double *xp, *fvecp, *fjac, *fvec, *x;
655 double *err;
656 int mode, m, n, ldfjac;
657
658 if (!PyArg_ParseTuple(args,"iiOOOiO!OiO!",&m, &n, &o_x, &o_fvec, &o_fjac, &ldfjac, &PyArray_Type, (PyObject **)&ap_xp, &o_fvecp, &mode, &PyArray_Type, (PyObject **)&ap_err)) return NULL;
659
660 ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(o_x,NPY_DOUBLE,1,1);
661 if (ap_x == NULL) goto fail;
662 if (n != PyArray_DIMS(ap_x)[0])
663 PYERR(minpack_error,"Input data array (x) must have length n");
664 x = (double *) PyArray_DATA(ap_x);
665 if (!PyArray_IS_C_CONTIGUOUS(ap_xp) || (PyArray_TYPE(ap_xp) != NPY_DOUBLE))
666 PYERR(minpack_error,"Seventh argument (xp) must be contiguous array of type Float64.");
667
668 if (mode == 1) {
669 fvec = NULL;
670 fjac = NULL;
671 xp = (double *)PyArray_DATA(ap_xp);
672 fvecp = NULL;
673 err = NULL;
674 CHKDER(&m, &n, x, fvec, fjac, &ldfjac, xp, fvecp, &mode, err);
675 }
676 else if (mode == 2) {
677 if (!PyArray_IS_C_CONTIGUOUS(ap_err) || (PyArray_TYPE(ap_err) != NPY_DOUBLE))
678 PYERR(minpack_error,"Last argument (err) must be contiguous array of type Float64.");
679 ap_fvec = (PyArrayObject *)PyArray_ContiguousFromObject(o_fvec,NPY_DOUBLE,1,1);
680 ap_fjac = (PyArrayObject *)PyArray_ContiguousFromObject(o_fjac,NPY_DOUBLE,2,2);
681 ap_fvecp = (PyArrayObject *)PyArray_ContiguousFromObject(o_fvecp,NPY_DOUBLE,1,1);
682 if (ap_fvec == NULL || ap_fjac == NULL || ap_fvecp == NULL) goto fail;
683
684 fvec = (double *)PyArray_DATA(ap_fvec);
685 fjac = (double *)PyArray_DATA(ap_fjac);
686 xp = (double *)PyArray_DATA(ap_xp);
687 fvecp = (double *)PyArray_DATA(ap_fvecp);
688 err = (double *)PyArray_DATA(ap_err);
689
690 CHKDER(&m, &n, x, fvec, fjac, &m, xp, fvecp, &mode, err);
691
692 Py_DECREF(ap_fvec);
693 Py_DECREF(ap_fjac);
694 Py_DECREF(ap_fvecp);
695 }
696 else
697 PYERR(minpack_error,"Invalid mode, must be 1 or 2.");
698
699 Py_DECREF(ap_x);
700
701 Py_INCREF(Py_None);
702 return Py_None;
703
704 fail:
705 Py_XDECREF(ap_fvec);
706 Py_XDECREF(ap_fjac);
707 Py_XDECREF(ap_fvecp);
708 Py_XDECREF(ap_x);
709 return NULL;
710 }
711
712
713
714
715
716
717
718
719
720