1 #include <pygsl/solver.h>
2 #include <pygsl/rng.h>
3 #include <gsl/gsl_monte.h>
4 #include <gsl/gsl_monte_plain.h>
5 #include <gsl/gsl_monte_miser.h>
6 #include <gsl/gsl_monte_vegas.h>
7
8 const char * filename = __FILE__;
9 PyObject *module = NULL;
10 static const char monte_f_type_name[] = "F-monte";
11 static const char module_doc[] = "XXX Missing";
12 static const char PyGSL_monte_type[] = "monte";
13
14 enum PyGSL_gsl_monte_type{
15 PyGSL_MONTE_plain = 1,
16 PyGSL_MONTE_miser = 2,
17 PyGSL_MONTE_vegas = 3
18 };
19
20 struct _monte_csys{
21 enum PyGSL_gsl_monte_type type;
22 };
23
24 typedef struct _monte_csys monte_csys;
25
26 static double
PyGSL_g_test(double * k,size_t dim,void * params)27 PyGSL_g_test(double *k, size_t dim, void *params)
28 {
29 double A = 1.0 / (M_PI * M_PI * M_PI);
30 return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
31 }
32
33 static double
PyGSL_monte_function_wrap(double * x,size_t dim,void * params)34 PyGSL_monte_function_wrap(double *x, size_t dim, void *params)
35 {
36 double tmp, *dresult;
37 int dimension;
38 PyGSL_solver *s;
39 PyArrayObject *a_array = NULL;
40 PyObject *result = NULL, *arglist=NULL, *callback, *retval;
41 PyGSL_error_info info;
42 /* the line number to appear in the traceback */
43 int trb_lineno = -1, i;
44 struct pygsl_array_cache * cache_ptr;
45 gsl_vector_view view;
46
47 FUNC_MESS_BEGIN();
48 s = (PyGSL_solver *) params;
49
50 /*
51 flag = PyGSL_function_wrap_On_O(&view.vector, s->cbs[0], s->args,
52 &tmp, NULL, view.vector.size,
53 "monte_wrap");
54 */
55
56 FUNC_MESS_BEGIN();
57 callback = s->cbs[0];
58 assert(s->args != NULL);
59 assert(callback != NULL);
60
61 /* Do I need to copy the array ??? */
62 for(i = 2; i < PyGSL_SOLVER_N_ARRAYS; ++i){
63 cache_ptr = &(s->cache[i]);
64 if(cache_ptr->ref == NULL)
65 /* no array found */
66 break;
67 if(x == cache_ptr->data){
68 a_array = cache_ptr->ref;
69 break;
70 }
71 }
72 if(i >= PyGSL_SOLVER_N_ARRAYS){
73 trb_lineno = __LINE__ -1;
74 pygsl_error("Not enough space to cache arrays...", filename, trb_lineno, GSL_ESANITY);
75 goto fail;
76 }
77 if(a_array == NULL){
78 dimension = dim;
79 a_array = (PyArrayObject *) PyArray_FromDimsAndData(1, &dimension, PyArray_DOUBLE, (char *)x);
80 cache_ptr->ref = a_array;
81 cache_ptr->data = x;
82 /*
83 * Required for deallocation.
84 * Tuples steal the references. But this array will be dereferenced by the
85 * normal procedure!
86 */
87 Py_INCREF(a_array);
88 }
89 if (a_array == NULL){
90 trb_lineno = __LINE__ - 2;
91 goto fail;
92 }
93 /* arglist was already set up */
94 result = (PyObject *) s->cache[1].ref;
95 dresult = (double *) s->cache[1].data;
96 arglist = (PyObject *) s->cache[0].ref;
97
98 Py_INCREF(s->args);
99 PyTuple_SET_ITEM(arglist, 0, (PyObject *) a_array);
100 PyTuple_SET_ITEM(arglist, 1, s->args);
101
102 /*
103 arglist = Py_BuildValue("(OOO)", a_array, s->args, result);
104 if(DEBUG > 2){
105 fprintf(stderr, "callback = %p, arglist = %p\n", callback, arglist);
106 }
107 */
108 FUNC_MESS(" Call Python Object BEGIN");
109 retval = PyEval_CallObject(callback, arglist);
110 FUNC_MESS(" Call Python Object END");
111
112 /*
113 info.callback = callback;
114 info.message = __FUNCTION__;
115 if(PyGSL_CHECK_PYTHON_RETURN(retval, 0, &info) != GSL_SUCCESS){
116 trb_lineno = __LINE__ - 1;
117 goto fail;
118 }
119 info.argnum = 1;
120 DEBUG_MESS(3, "result was %e", *dresult);
121 */
122 FUNC_MESS_END();
123 return *dresult;
124
125
126 fail:
127 PyGSL_add_traceback(NULL, __FILE__, __FUNCTION__, trb_lineno);
128 FUNC_MESS("Failure");
129
130 if(s->isset == 1){
131 longjmp(s->buffer, GSL_EFAILED);
132 FUNC_MESS("\t\t Using jump buffer");
133 } else {
134 FUNC_MESS("\t\t Jump buffer was not defined!");
135 }
136 tmp = gsl_nan();
137 return tmp;
138 }
139
140 static PyObject *
PyGSL_monte_init(PyGSL_solver * self,PyObject * args)141 PyGSL_monte_init(PyGSL_solver *self, PyObject *args)
142 {
143 int flag = GSL_EFAILED;
144 monte_csys * csys;
145
146 FUNC_MESS_BEGIN();
147 assert(PyGSL_solver_check(self));
148 csys = (monte_csys *)self->c_sys;
149 assert(csys);
150 switch(csys->type){
151 case PyGSL_MONTE_plain:
152 flag = gsl_monte_plain_init(self->solver);
153 break;
154 case PyGSL_MONTE_miser:
155 flag = gsl_monte_miser_init(self->solver);
156 break;
157 case PyGSL_MONTE_vegas:
158 flag = gsl_monte_vegas_init(self->solver);
159 break;
160 default:
161 DEBUG_MESS(2, "Monte type %d unknown",flag);
162 PyGSL_ERROR_NULL("Unknown monte type!", GSL_ESANITY);
163 }
164
165 if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
166 PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
167 return NULL;
168 }
169 Py_INCREF(Py_None);
170 FUNC_MESS_END();
171 return Py_None;
172 }
173
174 typedef int (*pygsl_monte_fptr_t)(gsl_monte_function * F, double XL[], double XU[], size_t DIM,
175 size_t CALLS, gsl_rng * R, gsl_monte_plain_state * S,
176 double * RESULT, double * ABSERR);
177
178 static PyObject *
PyGSL_monte_integrate(PyGSL_solver * self,PyObject * args)179 PyGSL_monte_integrate(PyGSL_solver *self, PyObject *args)
180 {
181
182 int flag = GSL_EFAILED, line, dim;
183 unsigned int ncalls;
184 monte_csys * csys;
185 gsl_rng * rng;
186 PyObject *func=NULL, *xlo=NULL, *xuo=NULL, *rng_obj=NULL, *ncallso=NULL, *argso=NULL;
187 PyArrayObject *xla=NULL;
188 PyArrayObject *xua=NULL;
189 double result, abserr;
190 gsl_monte_function mfunc;
191 pygsl_monte_fptr_t fptr;
192
193 FUNC_MESS_BEGIN();
194 assert(PyGSL_solver_check(self));
195
196 if(!PyArg_ParseTuple(args, "OOOOO|O", &func, &xlo, &xuo, &ncallso, &rng_obj, &argso)){
197 line = __LINE__ - 1;
198 goto fail;
199 }
200
201 if(!PyCallable_Check(func)){
202 line = __LINE__ - 1;
203 pygsl_error("The function argument must be callable!", filename, line, GSL_EINVAL);
204 goto fail;
205 }
206
207 if(PyGSL_PYLONG_TO_UINT(ncallso, &ncalls, NULL) != GSL_SUCCESS){
208 line = __LINE__ - 1;
209 goto fail;
210 }
211
212 if(!PyGSL_RNG_Check(rng_obj)){
213 line = __LINE__ - 1;
214 pygsl_error("The rng object must be a rng!", filename, line, GSL_EINVAL);
215 goto fail;
216 }
217
218 rng = ((PyGSL_rng *) rng_obj)->rng;
219
220 xla = PyGSL_vector_check(xlo, -1, PyGSL_DARRAY_CINPUT(2), NULL, NULL);
221 if(xla == NULL){
222 line = __LINE__ - 2;
223 goto fail;
224 }
225
226 dim = xla->dimensions[0];
227 xua = PyGSL_vector_check(xuo, dim, PyGSL_DARRAY_CINPUT(3), NULL, NULL);
228 if(xua == NULL){
229 line = __LINE__ - 2;
230 goto fail;
231 }
232
233 if(argso == NULL){
234 Py_INCREF(Py_None);
235 argso = Py_None;
236 }
237 Py_XDECREF(self->cbs[0]);
238 Py_INCREF(func);
239 self->cbs[0] = func;
240 self->args = argso;
241 csys = (monte_csys *)self->c_sys;
242
243 switch(csys->type){
244 case PyGSL_MONTE_plain: fptr = (pygsl_monte_fptr_t)&gsl_monte_plain_integrate; break;
245 case PyGSL_MONTE_miser: fptr = (pygsl_monte_fptr_t)&gsl_monte_miser_integrate; break;
246 case PyGSL_MONTE_vegas: fptr = (pygsl_monte_fptr_t)&gsl_monte_vegas_integrate; break;
247 default:
248 line = __LINE__ - 5;
249 DEBUG_MESS(2, "Monte type %d unknown",flag);
250 pygsl_error("Unknown monte type!", filename, line, GSL_ESANITY);
251 goto fail;
252 }
253
254 assert(fptr);
255 mfunc.f = &PyGSL_monte_function_wrap;
256 mfunc.f = &PyGSL_g_test;
257 mfunc.dim = dim;
258 mfunc.params = self;
259
260 if(setjmp(self->buffer) != 0){
261 self->isset = 0;
262 line = __LINE__ - 2;
263 goto fail;
264 }else{
265 self->isset = 1;
266 }
267 flag = fptr(&mfunc, (double *) xla->data, (double *)xua->data,
268 dim, ncalls, rng, self->solver, &result, &abserr);
269 self->isset = 0;
270 if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
271 line = __LINE__ - 2;
272 return NULL;
273 }
274
275 Py_DECREF(xla); xla = NULL;
276 Py_DECREF(xua); xua = NULL;
277 Py_DECREF(self->cbs[0]);
278 Py_DECREF(self->args);
279 self->cbs[0] = NULL;
280 self->args = NULL;
281 FUNC_MESS_END();
282
283 return Py_BuildValue("dd", result, abserr);
284
285 fail:
286 FUNC_MESS("FAIL");
287 Py_XDECREF(xla);
288 Py_XDECREF(xua);
289 Py_XDECREF(self->cbs[0]);
290 Py_XDECREF(self->args);
291 self->cbs[0] = NULL;
292 self->args = NULL;
293 PyGSL_add_traceback(module, filename, __FUNCTION__, line);
294 return NULL;
295 }
296
297 #define GET_SET(montetype, name, mode) \
298 static PyObject * GetSet_ ## montetype ## _ ## name(PyGSL_solver *self, PyObject *args) \
299 { \
300 gsl_monte_## montetype ## _state *state = self->solver; \
301 return PyGSL_solver_GetSet((PyObject *)self, args, &(state->name), mode); \
302 }
303
304 #define GET_SET_DOUBLE(mintype, name) GET_SET(mintype, name, PyGSL_MODE_DOUBLE)
305 #define GET_SET_SIZE_T(mintype, name) GET_SET(mintype, name, PyGSL_MODE_SIZE_T)
306 #define GET_SET_INT(mintype, name) GET_SET(mintype, name, PyGSL_MODE_INT)
307
308 GET_SET_DOUBLE(miser, estimate_frac)
309 GET_SET_SIZE_T(miser, min_calls)
310 GET_SET_SIZE_T(miser, min_calls_per_bisection)
311 GET_SET_DOUBLE(miser, alpha)
312 GET_SET_DOUBLE(miser, dither)
313
314 GET_SET_DOUBLE(vegas, result)
315 GET_SET_DOUBLE(vegas, sigma)
316 GET_SET_DOUBLE(vegas, chisq)
317 GET_SET_DOUBLE(vegas, alpha)
318 GET_SET_SIZE_T(vegas, iterations)
319 GET_SET_INT(vegas, stage)
320 GET_SET_INT(vegas, mode)
321 GET_SET_INT(vegas, verbose)
322
323 #undef GET_SET_DOUBLE
324 #undef GET_SET_SIZE_T
325 #undef GET_SET_INT
326 #undef GET_SET
327 #define GETSET(montetype, name, mode) {#name, (PyCFunction) GetSet_ ## montetype ## _ ## name, METH_VARARGS, NULL},
328 #define GET_SET_DOUBLE(montetype, name) GETSET(montetype, name, 0)
329 #define GET_SET_INT(montetype, name) GETSET(montetype, name, 0)
330 #define GET_SET_SIZE_T(montetype, name) GETSET(montetype, name, 0)
331
332 #define MONTE_STANDARD_METHODS \
333 {"init", (PyCFunction)PyGSL_monte_init, METH_VARARGS, NULL}, \
334 {"integrate", (PyCFunction)PyGSL_monte_integrate, METH_VARARGS, NULL}, \
335
336 static PyMethodDef PyGSL_monte_plain_methods[] = {
337 MONTE_STANDARD_METHODS
338 {NULL, NULL, 0, NULL}
339 };
340
341 static PyMethodDef PyGSL_monte_miser_methods[] = {
342 MONTE_STANDARD_METHODS
343 GET_SET_DOUBLE(miser, estimate_frac)
344 GET_SET_SIZE_T(miser, min_calls)
345 GET_SET_SIZE_T(miser, min_calls_per_bisection)
346 GET_SET_DOUBLE(miser, alpha)
347 GET_SET_DOUBLE(miser, dither)
348 {NULL, NULL, 0, NULL}
349 };
350
351 static PyMethodDef PyGSL_monte_vegas_methods[] = {
352 MONTE_STANDARD_METHODS
353 GET_SET_DOUBLE(vegas, result)
354 GET_SET_DOUBLE(vegas, sigma)
355 GET_SET_DOUBLE(vegas, chisq)
356 GET_SET_DOUBLE(vegas, alpha)
357 GET_SET_SIZE_T(vegas, iterations)
358 GET_SET_INT(vegas, stage)
359 GET_SET_INT(vegas, mode)
360 GET_SET_INT(vegas, verbose)
361 {NULL, NULL, 0, NULL}
362 };
363
364
365 /* make the init look similar to the one */
366 static void *
PyGSL_gsl_monte_alloc(void * type,int n)367 PyGSL_gsl_monte_alloc(void * type, int n)
368 {
369 enum PyGSL_gsl_monte_type flag = (enum PyGSL_gsl_monte_type) type;
370 void * result = NULL;
371
372 FUNC_MESS_BEGIN();
373 switch(flag){
374 case PyGSL_MONTE_plain:
375 result = gsl_monte_plain_alloc(n);
376 break;
377 case PyGSL_MONTE_miser:
378 result = gsl_monte_miser_alloc(n);
379 break;
380 case PyGSL_MONTE_vegas:
381 result = gsl_monte_vegas_alloc(n);
382 break;
383 default:
384 DEBUG_MESS(2, "Monte type %d unknown",flag);
385 PyGSL_ERROR_NULL("Unknown monte type!", GSL_ESANITY);
386 }
387 FUNC_MESS_END();
388 return result;
389 }
390
391 const struct _SolverStatic
392 monte_plain_solver_f = {{(void_m_t) gsl_monte_plain_free,
393 (void_m_t) NULL,
394 NULL,
395 NULL},
396 1, PyGSL_monte_plain_methods, PyGSL_monte_type},
397 monte_miser_solver_f = {{(void_m_t) gsl_monte_miser_free,
398 (void_m_t) NULL,
399 NULL,
400 NULL},
401 1, PyGSL_monte_miser_methods, PyGSL_monte_type},
402 monte_vegas_solver_f = {{(void_m_t) gsl_monte_vegas_free,
403 (void_m_t) NULL,
404 NULL,
405 NULL},
406 1, PyGSL_monte_vegas_methods, PyGSL_monte_type};
407
408 static PyObject*
PyGSL_init_monte(PyObject * self,PyObject * args,int type,const struct _SolverStatic * sost)409 PyGSL_init_monte(PyObject *self, PyObject *args, int type, const struct _SolverStatic *sost)
410 {
411 PyGSL_solver *result = NULL;
412 struct pygsl_array_cache *tmp;
413 int line = -1, dims;
414 monte_csys * csys;
415 PyObject *mytuple;
416
417 FUNC_MESS_BEGIN();
418 const solver_alloc_struct alloc = {(const void *) type, (void *) PyGSL_gsl_monte_alloc, sost};
419 result = (PyGSL_solver *) PyGSL_solver_dn_init(self, args, &alloc, 1);
420 if(result == NULL){
421 line = __LINE__ - 2;
422 goto fail;
423 }
424
425 csys = (monte_csys *) calloc(1, sizeof(monte_csys));
426 if(csys == NULL){
427 PyErr_NoMemory();
428 line = __LINE__ - 1;
429 goto fail;
430 }
431
432
433 dims = 1;
434 tmp = result->cache;
435
436 tmp[1].ref = PyGSL_New_Array(1, &dims, PyArray_DOUBLE);
437 tmp[1].data = (double *) tmp[1].ref->data;
438 mytuple = PyTuple_New(3);
439 PyTuple_SetItem(mytuple, 2, (PyObject *) tmp[1].ref);
440 /*
441 * Required for deallocation.
442 * Tuples steal the references. But this array will be dereferenced by the
443 * normal procedure!
444 */
445 Py_INCREF(tmp[1].ref);
446 Py_INCREF(tmp[1].ref);
447
448 /*
449 * Initalise the others
450 */
451 Py_INCREF(Py_None);
452 Py_INCREF(Py_None);
453 Py_INCREF(Py_None);
454 Py_INCREF(Py_None);
455 PyTuple_SetItem(mytuple, 0, Py_None);
456 PyTuple_SetItem(mytuple, 1, Py_None);
457
458 tmp[0].ref = (PyArrayObject*) mytuple;
459
460 result->cache = tmp;
461 csys->type = type;
462 result->c_sys = csys;
463
464 FUNC_MESS_END();
465 return (PyObject *) result;
466
467 fail:
468 FUNC_MESS("Fail");
469 Py_XDECREF(result);
470 PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__);
471 return NULL;
472 }
473
474 #define MONTE_INIT(name) \
475 static PyObject* PyGSL_init_ ## name(PyObject * self, PyObject *args) \
476 { return PyGSL_init_monte(self, args, PyGSL_MONTE_ ##name, &monte_ ## name ## _solver_f); }
477
478 MONTE_INIT(plain)
479 MONTE_INIT(miser)
480 MONTE_INIT(vegas)
481
482
483 static PyMethodDef mMethods[] = {
484 {"plain", PyGSL_init_plain, METH_VARARGS, NULL},
485 {"miser", PyGSL_init_miser, METH_VARARGS, NULL},
486 {"vegas", PyGSL_init_vegas, METH_VARARGS, NULL},
487 {NULL, NULL, 0, NULL}
488 };
489
490 void
initmonte(void)491 initmonte(void)
492 {
493 PyObject* m, *dict, *item;
494 FUNC_MESS_BEGIN();
495
496 m=Py_InitModule("monte", mMethods);
497 module = m;
498 assert(m);
499 dict = PyModule_GetDict(m);
500 if(!dict)
501 goto fail;
502
503 init_pygsl()
504 import_pygsl_solver();
505 assert(PyGSL_API);
506
507
508 if (!(item = PyString_FromString((char*)module_doc))){
509 PyErr_SetString(PyExc_ImportError,
510 "I could not generate module doc string!");
511 goto fail;
512 }
513 if (PyDict_SetItemString(dict, "__doc__", item) != 0){
514 PyErr_SetString(PyExc_ImportError,
515 "I could not init doc string!");
516 goto fail;
517 }
518 PyModule_AddIntConstant(m, "GSL_VEGAS_MODE_IMPORTANCE", GSL_VEGAS_MODE_IMPORTANCE);
519 PyModule_AddIntConstant(m, "GSL_VEGAS_MODE_IMPORTANCE_ONLY", GSL_VEGAS_MODE_IMPORTANCE_ONLY);
520 PyModule_AddIntConstant(m, "GSL_VEGAS_MODE_STRATIFIED", GSL_VEGAS_MODE_STRATIFIED);
521
522 FUNC_MESS_END();
523 return;
524
525 fail:
526 FUNC_MESS("FAIL");
527 return;
528
529 }
530