1 /* -*- C -*- */
2 /**
3 * author: Pierre Schnizer
4 */
5
6 /*
7 * The gsl_multifit_linear_workspace structre was changed from GSL 1.X to
8 * GSL 2.X.
9 * setup checks for that and then appropriate measures are taken
10 */
11 %{
12 #include <gsl/gsl_multifit.h>
13 #include <gsl/gsl_fit.h>
14 #include <gsl/gsl_matrix.h>
15 #include <gsl/gsl_vector.h>
16 #include <pygsl/pygsl_features.h>
17
18 #define _PyGSL_TO_ARRAY_INDEX_CAST(arg) ((PyGSL_array_index_t) (arg))
19 #if PYGSL_GSL_MAJOR_VERSION == 2
20 #if PYGSL_GSL_MINOR_VERSION == 0
21 #error GSL 2.0 implementation of multifit has a bug in multifit_wliner
22 #endif /* PYGSL_GSL_MINOR_VERSION == 0 */
23 #endif /* PYGSL_GSL_MAJOR_VERSION == 2 */
24
25 #ifdef _PYGSL_GSL_HAS_MULTFIT_LINEAR_WORKSPACE_STRUCT_MEMBER_NMAX_PMAX
26 #define _PyGSL_MULTIFIT_GET_NMAX(s) _PyGSL_TO_ARRAY_INDEX_CAST((s->nmax))
27 #define _PyGSL_MULTIFIT_GET_PMAX(s) _PyGSL_TO_ARRAY_INDEX_CAST((s->pmax))
28 /* fix a bug in multifit_wlinear ... */
29 #define _PyGSL_MULTIFIT_WORKSPACE_SET(s) s->n = s->nmax; (s->p = s->pmax);
30 #else /* _PYGSL_GSL_HAS_MULTFIT_LINEAR_WORKSPACE_STRUCT_MEMBER_NMAX_PMAX */
31 #ifndef _PYGSL_GSL_HAS_MULTFIT_LINEAR_WORKSPACE_STRUCT_MEMBER_N_P
32 #error "I expected that if nmax and pmax is not found that n and p are always"
33 #error "found"
34 #endif /* _PYGSL_GSL_HAS_MULTFIT_LINEAR_WORKSPACE_STRUCT_MEMBER_N_P */
35 #define _PyGSL_MULTIFIT_WORKSPACE_SET(s)
36 #define _PyGSL_MULTIFIT_GET_NMAX(s) _PyGSL_TO_ARRAY_INDEX_CAST((s->n))
37 #define _PyGSL_MULTIFIT_GET_PMAX(s) _PyGSL_TO_ARRAY_INDEX_CAST((s->p))
38 #endif /* _PYGSL_GSL_HAS_MULTFIT_LINEAR_WORKSPACE_STRUCT_MEMBER_NMAX_PMAX */
39 %}
40 %include typemaps.i
41 %include gsl_block_typemaps.i
42
43 %typemap(arginit) const gsl_matrix * IN_AND_SIZE %{
44 PyArrayObject * _PyMatrix$argnum = NULL;
45 TYPE_VIEW_$1_basetype _matrix$argnum;
46 PyGSL_array_index_t _mat_dim0_arg$argnum = -1;
47 PyGSL_array_index_t _mat_dim1_arg$argnum = -1;
48
49 %}
50 %typemap(check) const gsl_matrix * IN_AND_SIZE %{
51 assert(_PyMatrix$argnum != NULL);
52 _mat_dim0_arg$argnum = PyArray_DIM(_PyMatrix$argnum, 0);
53 _mat_dim1_arg$argnum = PyArray_DIM(_PyMatrix$argnum, 1);
54 %}
55
56 %typemap(in) gsl_matrix * IN_AND_SIZE = gsl_matrix * IN;
57 %typemap(freearg) gsl_matrix * IN_AND_SIZE = gsl_matrix * IN;
58
59
60 gsl_multifit_linear_workspace *
61 gsl_multifit_linear_alloc (size_t n, size_t p);
62
63 void
64 gsl_multifit_linear_free (gsl_multifit_linear_workspace * work);
65
66 %typemap(arginit) gsl_multifit_linear_workspace * work_provide %{
67 PyGSL_array_index_t _work_provide_n_$1_name = -1;
68 PyGSL_array_index_t _work_provide_p_$1_name = -1;
69 %}
70
71 %typemap( in) gsl_multifit_linear_workspace * work_provide {
72 if ((SWIG_ConvertPtr($input, (void **) &$1, SWIGTYPE_p_gsl_multifit_linear_workspace,1)) == -1){
73 goto fail;
74 }
75 _work_provide_n_$1_name = _PyGSL_MULTIFIT_GET_NMAX($1);
76 _work_provide_p_$1_name = _PyGSL_MULTIFIT_GET_PMAX($1);
77 _PyGSL_MULTIFIT_WORKSPACE_SET($1);
78 DEBUG_MESS(2, "work->n = %ld work ->p = %ld",
79 (long) _work_provide_n_$1_name, (long) _work_provide_p_$1_name);
80 /*
81 DEBUG_MESS(2, "work->n = %ld work ->p = %ld work->nmax = %ld, work->pmax = %ld",
82 (long) _work_provide_n_$1_name, (long) _work_provide_p_$1_name,
83 (long) $1->nmax, (long) $1->pmax);
84 */
85 };
86
87 %typemap( in, numinputs = 0) gsl_vector * OUT %{
88 /* All done in check as the workspace stores the information about the required size */
89 %}
90
91 %typemap(check) gsl_vector * OUT {
92 PyGSL_array_index_t stride, lvec;
93 lvec = _work_provide_p_work_provide;
94 lvec = _mat_dim1_arg1;
95
96 DEBUG_MESS(2, "out vector length = %ld not %ld", (long) lvec, (long) _work_provide_p_work_provide);
97 _PyVector$argnum = (PyArrayObject *) PyGSL_New_Array(1, &lvec, NPY_DOUBLE);
98 if(NULL == _PyVector$argnum){
99 goto fail;
100 }
101
102 if(PyGSL_STRIDE_RECALC(PyArray_STRIDE(_PyVector$argnum, 0), sizeof(BASIS_TYPE($1_basetype)), &stride) != GSL_SUCCESS)
103 goto fail;
104
105 _vector$argnum = TYPE_VIEW_ARRAY_STRIDES_$1_basetype((BASIS_C_TYPE($1_basetype) *) PyArray_DATA(_PyVector$argnum),
106 stride,
107 PyArray_DIM(_PyVector$argnum, 0));
108 $1 = ($basetype *) &(_vector$argnum.vector);
109
110 }
111 %typemap( in) gsl_matrix * OUT = gsl_vector * OUT;
112 %typemap(check) gsl_matrix * OUT {
113 PyArrayObject * a_array;
114 BASIS_C_TYPE($1_basetype) *data = NULL;
115 PyGSL_array_index_t stride_recalc=0, dimensions[2];
116
117
118 dimensions[0] = _mat_dim1_arg1;
119 dimensions[1] = _mat_dim1_arg1;
120 a_array = (PyArrayObject *) PyGSL_New_Array(2, dimensions, NPY_DOUBLE);
121 if(NULL == a_array){
122 goto fail;
123 }
124 _PyMatrix$argnum = a_array;
125
126
127 if(PyGSL_STRIDE_RECALC(PyArray_STRIDE(a_array, 0), sizeof(BASIS_TYPE($1_basetype)), &stride_recalc) != GSL_SUCCESS)
128 goto fail;
129 /* (BASIS_TYPE_$1_basetype *) */
130 data = (BASIS_C_TYPE($1_basetype) *) PyArray_DATA(a_array);
131 assert(data != NULL);
132 _matrix$argnum = TYPE_VIEW_ARRAY_$1_basetype(data, PyArray_DIM(a_array,0), PyArray_DIM(a_array, 1));
133 assert(_matrix$argnum.matrix.data != NULL);
134 $1 = ($basetype *) &(_matrix$argnum.matrix);
135 DEBUG_MESS(2, "matrix: data %p size = [%ld, %ld]", (void *) _matrix$argnum.matrix.data,
136 (long) _matrix$argnum.matrix.size1,
137 (long) _matrix$argnum.matrix.size2
138 );
139
140 }
141 %typemap(argout) gsl_vector * OUT{
142 $result = SWIG_Python_AppendOutput($result, (PyObject *) _PyVector$argnum);
143 _PyVector$argnum =NULL;
144 };
145 %typemap(argout) gsl_matrix * OUT{
146 $result = SWIG_Python_AppendOutput($result, (PyObject *) _PyMatrix$argnum);
147 _PyMatrix$argnum =NULL;
148 };
149
150 /* ---------------------------------------------------------------------- */
151 /*
152 * A typemap for translating numpy arrays to the least squares
153 * representation.
154 *
155 * The length of vector x establishes the reference. All others check
156 * their length to this vector. So if you apply this typemap it will
157 * fail if no variable was named x!!!
158 *
159 *
160 */
161
162 %typemap(arginit) (double *, size_t ) %{
163 PyArrayObject *_PyVector$argnum = NULL;
164 size_t _PyVectorLength$1_name = 0;
165 %}
166
167 %typemap(in) (double * , size_t ){
168 PyGSL_array_index_t strides, size;
169 /* This should be a preprocessor directive. */
170 if ( '$1_name' == 'x' )
171 size = -1;
172 else
173 size = _PyVectorLengthx;
174 _PyVector$argnum = PyGSL_vector_check($input, size, PyGSL_DARRAY_INPUT($argnum), &strides, NULL);
175 if (_PyVector$argnum == NULL)
176 goto fail;
177
178 $1 = (double *) PyArray_DATA(_PyVector$argnum);
179 $2 = (size_t) strides;
180 _PyVectorLength$1_name = (size_t) PyArray_DIM(_PyVector$argnum, 0);
181 };
182 %typemap(argout) (double *, size_t ) {
183 Py_XDECREF(_PyVector$argnum);
184 };
185
186 %typemap( in, numinputs=0) const size_t n {
187 $1 = 0;
188 };
189
190 %typemap(check) const size_t n {
191 $1 = _PyVectorLengthx;
192 };
193
194 %apply double * OUTPUT {double * c0,
195 double * c1,
196 double * cov00,
197 double * cov01,
198 double * cov11,
199 double * sumsq,
200 double * chisq,
201 double * Y,
202 double * Y_ERR};
203
204 gsl_error_flag_drop
205 gsl_multifit_linear (const gsl_matrix * IN_AND_SIZE,
206 const gsl_vector * IN,
207 gsl_vector * OUT,
208 gsl_matrix * OUT,
209 double * chisq,
210 gsl_multifit_linear_workspace * work_provide);
211
212 #if 0
213 #if PyGSL_GSL_MAJOR_VERSION >= 2
214 /* needs works space ... implement ... */
215 #error "Work space needs to be correctly implemented "
216 gsl_error_flag_drop
217 gsl_multifit_linear_svd (const gsl_matrix * IN,
218 gsl_multifit_linear_workspace * work_provide);
219 gsl_error_flag_drop
220 gsl_multifit_linear_bsvd (const gsl_matrix * IN,
221 gsl_multifit_linear_workspace * work_provide);
222 #else /* PyGSL_GSL_MAJOR_VERSION >= 2 */
223 gsl_error_flag_drop
224 gsl_multifit_linear_svd (const gsl_matrix * IN,
225 const gsl_vector * IN,
226 double TOL,
227 size_t * OUTPUT,
228 gsl_vector * OUT,
229 gsl_matrix * OUT2,
230 double * chisq,
231 gsl_multifit_linear_workspace * work_provide);
232 #endif /* PyGSL_GSL_MAJOR_VERSION >= 2 */
233 #endif
234
235 %apply gsl_vector *IN {gsl_vector *y};
236 %apply gsl_vector *IN {gsl_vector *w};
237 %apply gsl_vector *IN {gsl_vector *c};
238 %apply gsl_vector *IN {gsl_vector *r};
239 %apply gsl_matrix *IN {gsl_matrix *x};
240 %apply gsl_matrix *IN {gsl_matrix *cov};
241 gsl_error_flag_drop
242 gsl_multifit_wlinear (const gsl_matrix * IN_AND_SIZE,
243 const gsl_vector * w,
244 const gsl_vector * y,
245 gsl_vector * OUT,
246 gsl_matrix * OUT,
247 double * chisq,
248 gsl_multifit_linear_workspace * work_provide);
249
250
251 gsl_error_flag_drop
252 gsl_multifit_wlinear_svd (const gsl_matrix * IN_AND_SIZE,
253 const gsl_vector * w,
254 const gsl_vector * y,
255 double TOL,
256 size_t * OUTPUT,
257 gsl_vector * OUT,
258 gsl_matrix * OUT,
259 double * chisq,
260 gsl_multifit_linear_workspace * work_provide);
261
262 gsl_error_flag_drop
263 gsl_multifit_linear_est (const gsl_vector * x,
264 const gsl_vector * c,
265 const gsl_matrix * IN,
266 double * Y, double *Y_ERR);
267
268 /*
269 gsl_error_flag_drop
270 gsl_multifit_linear_residuals(const gsl_matrix * IN,
271 const gsl_vector * y,
272 const gsl_vector * c,
273 gsl_vector * r);
274
275 */
276
277
278 PyObject *
279 gsl_multifit_linear_est_matrix (const gsl_matrix * x,
280 const gsl_vector * c,
281 const gsl_matrix * cov);
282
283 PyObject *
284 pygsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
285 const gsl_vector *c);
286
287 %{
288 PyObject *
pygsl_multifit_linear_residuals(const gsl_matrix * X,const gsl_vector * y,const gsl_vector * c)289 pygsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
290 const gsl_vector *c)
291 {
292
293 int flag, line = __LINE__;
294
295 PyArrayObject *r_a = NULL;
296 gsl_vector_view r;
297 PyGSL_array_index_t dim = 0;
298 FUNC_MESS_BEGIN();
299
300 dim = y->size;
301 r_a = PyGSL_New_Array(1, &dim, NPY_DOUBLE);
302 if(r_a == NULL){
303 goto fail;
304 }
305 r = gsl_vector_view_array(PyArray_DATA(r_a), PyArray_DIM(r_a, 0));
306 flag = gsl_multifit_linear_residuals(X, y, c, &r.vector);
307 if(GSL_SUCCESS != PyGSL_ERROR_FLAG(flag)){
308 goto fail;
309 }
310 FUNC_MESS_END();
311
312 return (PyObject *) r_a;
313 fail:
314 FUNC_MESS("Fail");
315 Py_XDECREF(r_a);
316 return NULL;
317 }
318 %}
319
320 %apply (double *, size_t){(const double * x, const size_t xstride),
321 (const double * y, const size_t ystride),
322 (const double * w, const size_t wstride)
323 };
324 /* ---------------------------------------------------------------------- */
325 gsl_error_flag_drop gsl_fit_linear (const double * x, const size_t xstride,
326 const double * y, const size_t ystride,
327 const size_t n,
328 double * c0, double * c1,
329 double * cov00, double * cov01, double * cov11,
330 double * sumsq);
331
332
333 gsl_error_flag_drop gsl_fit_wlinear (const double * x, const size_t xstride,
334 const double * w, const size_t wstride,
335 const double * y, const size_t ystride,
336 const size_t n,
337 double * c0, double * c1,
338 double * cov00, double * cov01, double * cov11,
339 double * chisq);
340
341 %apply double * OUTPUT {double * y,
342 double * y_err}
343
344 gsl_error_flag_drop
345 gsl_fit_linear_est (const double x,
346 const double c0, const double c1,
347 const double c00, const double c01, const double c11,
348 double *y, double *y_err);
349
350
351 gsl_error_flag_drop gsl_fit_mul (const double * x, const size_t xstride,
352 const double * y, const size_t ystride,
353 const size_t n,
354 double * c1,
355 double * cov11,
356 double * sumsq);
357
358 gsl_error_flag_drop gsl_fit_wmul (const double * x, const size_t xstride,
359 const double * w, const size_t wstride,
360 const double * y, const size_t ystride,
361 const size_t n,
362 double * c1,
363 double * cov11,
364 double * sumsq);
365
366
367 gsl_error_flag_drop
368 gsl_fit_mul_est (const double x,
369 const double c1,
370 const double c11,
371 double *y, double *y_err);
372
373
374 /*
375 int gsl_fit_poly (const double * x,
376 const double * w,
377 const double * y,
378 size_t n,
379 double * c, size_t m,
380 double * chisq);
381
382 int gsl_fit_fns (const double * A,
383 const double * w,
384 const double * y,
385 size_t n,
386 double * c, size_t m,
387 double * chisq);
388
389 int gsl_fit_linear_nd (double * m, double * y, double * w);
390 */
391