← Back to branch summary

~nipy-developers/nipy/fff2

« back to all changes in this revision

Viewing changes to python/misc/routines.c

  • Committer: Bertrand THIRION
  • Date: 2008-04-03 17:29:55 UTC
  • mfrom: (13.1.5 fff2)
  • Revision ID: bt206016@is143015-20080403172955-w7v1pdjdriofvzzs
Merged Alexis'changes

Show diffs side-by-side

added added

removed removed

Lines of Context:
python/misc/routines.c
 
1
#include <fffpy.h>
 
2
#include <fff.h>
 
3
#include <fff_vector.h>
 
4
#include <fff_gen_stats.h>
 
5
 
 
6
 
 
7
#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API
 
8
 
 
9
#define DEFAULT_AXIS 0
 
10
 
 
11
/* Code pour la creation du module */ 
 
12
static char quantile_doc[] = 
 
13
" q = quantile(data, ratio=.5, interp=False, axis=0)\n\
 
14
Comments to follow";
 
15
static char mahalanobis_doc[] = 
 
16
" d2 = mahalanobis(data, vardata)\n\
 
17
data and vardata assumed C contiguous. More comments to follow";
 
18
static char module_doc[] = 
 
19
" Basic statistical routines.\n\
 
20
Author: Alexis Roche (CEA, Orsay, France), 2004-2007.";
 
21
 
 
22
 
 
23
/* 
 
24
   Quantile
 
25
 
 
26
*/
 
27
 
 
28
#define DEFAULT_RATIO 0.5
 
29
#define DEFAULT_INTERP 0 
 
30
 
 
31
static PyObject* quantile(PyObject* self, PyObject* args, PyObject* kw) 
 
32
{
 
33
  double ratio = DEFAULT_RATIO;
 
34
  int interp = DEFAULT_INTERP;
 
35
  int axis = DEFAULT_AXIS;
 
36
  PyArrayObject *X, *Y; 
 
37
  fff_vector *x, *y; 
 
38
  PyArrayMultiIterObject* multi;
 
39
  npy_intp ndim, dim_axis;
 
40
  npy_intp* dimAux; 
 
41
  static char *kwlist[] = {"data", "ratio", "interp", "axis", NULL}; 
 
42
 
 
43
  /* Parse input */ 
 
44
  if ( !PyArg_ParseTupleAndKeywords( args, kw, "O!|dii:quantile", kwlist, 
 
45
                                     &PyArray_Type, &X, &ratio, &interp, &axis ) )
 
46
    Py_RETURN_NONE; 
 
47
  
 
48
  /* Get input dimensions */
 
49
  ndim = PyArray_NDIM( X ); 
 
50
  dim_axis = PyArray_DIM( X, axis );
 
51
  dimAux = (npy_intp*) malloc( ndim*sizeof(npy_intp) );
 
52
  memcpy( (void*)dimAux, (void*)PyArray_DIMS(X), ndim*sizeof(npy_intp) ); 
 
53
 
 
54
  /* Allocate output array */
 
55
  dimAux[axis] = 1;
 
56
  Y = (PyArrayObject*) PyArray_SimpleNew(ndim, dimAux, NPY_DOUBLE); 
 
57
 
 
58
  /* Create multi iterator */ 
 
59
  multi = (PyArrayMultiIterObject*)PyArray_MultiIterAllButAxis(2, axis, X, Y);
 
60
 
 
61
  /* Create fff vectors (copy or not) */ 
 
62
  x = fff_vector_fromPyArrayIter(multi->iters[0], axis); 
 
63
  y = fff_vector_fromPyArrayIter(multi->iters[1], axis); 
 
64
 
 
65
  /* DEBUG */
 
66
  fprintf(stderr, "multi->index=%d  multi->size=%d\n", multi->index, multi->size); 
 
67
 
 
68
  /* Loop */ 
 
69
  while(multi->index < multi->size) {
 
70
 
 
71
    /* DEBUG */
 
72
    fprintf(stderr, "Everything fine.\n"); 
 
73
 
 
74
    fff_vector_get_fromPyArrayIter(x, multi->iters[0], axis); 
 
75
    fff_vector_get_fromPyArrayIter(y, multi->iters[1], axis); 
 
76
 
 
77
    /* DEBUG */
 
78
    fprintf(stderr, "Everything okay.\n"); 
 
79
  
 
80
    y->data[0] = fff_vector_quantile(x, ratio, interp); 
 
81
 
 
82
    /* DEBUG */
 
83
    fprintf(stderr, "So fine...\n"); 
 
84
 
 
85
    fff_vector_const_toPyArrayIter(y, multi->iters[1], axis); 
 
86
 
 
87
    /* DEBUG */
 
88
    fprintf(stderr, "Hey hey... still alive!\n"); 
 
89
 
 
90
    PyArray_MultiIter_NEXT(multi);
 
91
 
 
92
  }
 
93
 
 
94
  /* Free memory */ 
 
95
  free(dimAux); 
 
96
  fff_vector_delete(x);
 
97
  fff_vector_delete(y);
 
98
  Py_XDECREF(multi);
 
99
  
 
100
  /* Output */
 
101
  return (PyObject*)Y;
 
102
}
 
103
 
 
104
 
 
105
 
 
106
static PyObject* mahalanobis(PyObject* self, PyObject* args, PyObject* kw) 
 
107
{
 
108
  PyArrayObject *X, *VX, *D2; 
 
109
  fff_vector *x, *vx, *d2; 
 
110
  PyArrayMultiIterObject* multi;
 
111
  int axis = DEFAULT_AXIS;
 
112
  npy_intp ndim, n;
 
113
  npy_intp* dimAux;  
 
114
  fff_vector *x_tmp, *vx_tmp;
 
115
  fff_matrix *Sf;  
 
116
  fff_matrix Sview; 
 
117
  static char *kwlist[] = {"data", "vardata", NULL}; 
 
118
 
 
119
  /* Parse input */ 
 
120
  if ( !PyArg_ParseTupleAndKeywords( args, kw, "O!O!:mahalanobis", kwlist, 
 
121
                                     &PyArray_Type, &X,
 
122
                                     &PyArray_Type, &VX ) )
 
123
    Py_RETURN_NONE; 
 
124
  
 
125
  
 
126
  /* Dimensions */
 
127
  n = PyArray_DIM(X, 0); 
 
128
  ndim = PyArray_NDIM(X); 
 
129
  
 
130
  /* Allocate output array */
 
131
  dimAux = (npy_intp*) malloc( ndim*sizeof(npy_intp) );
 
132
  memcpy( (void*)dimAux, (void*)PyArray_DIMS(X), ndim*sizeof(npy_intp) ); 
 
133
  dimAux[0] = 1;
 
134
  D2 = (PyArrayObject*) PyArray_SimpleNew( ndim, dimAux, NPY_DOUBLE ); 
 
135
  
 
136
  /* Allocate local structures */ 
 
137
  x_tmp = fff_vector_new(n); 
 
138
  vx_tmp = fff_vector_new(n*n); 
 
139
  Sf = fff_matrix_new(n, n); 
 
140
  
 
141
  /* Create multi iterator */ 
 
142
  multi = (PyArrayMultiIterObject*)PyArray_MultiIterAllButAxis(3, axis, X, VX, D2);
 
143
  
 
144
  /* Create fff vectors (copy or not) */ 
 
145
  x = fff_vector_fromPyArrayIter(multi->iters[0], axis); 
 
146
  vx = fff_vector_fromPyArrayIter(multi->iters[1], axis); 
 
147
  d2 = fff_vector_fromPyArrayIter(multi->iters[2], axis); 
 
148
  
 
149
  /* Loop */ 
 
150
  while(multi->index < multi->size) {
 
151
 
 
152
    fff_vector_get_fromPyArrayIter(x, multi->iters[0], axis); 
 
153
    fff_vector_get_fromPyArrayIter(vx, multi->iters[1], axis); 
 
154
    fff_vector_get_fromPyArrayIter(d2, multi->iters[2], axis); 
 
155
 
 
156
    fff_vector_memcpy(x_tmp, x); 
 
157
    fff_vector_memcpy(vx_tmp, vx); 
 
158
    Sview = fff_matrix_view(vx_tmp->data, n, n, n); /* OK because vx_tmp is contiguous */ 
 
159
    d2->data[0] = fff_mahalanobis(x_tmp, &Sview, Sf);
 
160
    
 
161
    fff_vector_const_toPyArrayIter(d2, multi->iters[2], axis); 
 
162
    PyArray_MultiIter_NEXT(multi);
 
163
 
 
164
  }
 
165
 
 
166
  /* Free memory */ 
 
167
  free( dimAux ); 
 
168
  fff_vector_delete(x_tmp); 
 
169
  fff_vector_delete(vx_tmp); 
 
170
  fff_matrix_delete(Sf); 
 
171
  Py_XDECREF(multi);
 
172
  
 
173
  /* Output */
 
174
  return (PyObject*)D2; 
 
175
}
 
176
 
 
177
 
 
178
 
 
179
 
 
180
 
 
181
static PyMethodDef module_methods[] = {
 
182
  {"quantile", (PyCFunction)quantile, METH_KEYWORDS, quantile_doc},
 
183
  {"mahalanobis", (PyCFunction)mahalanobis, METH_KEYWORDS, mahalanobis_doc},
 
184
  {NULL, NULL, 0, NULL}
 
185
};
 
186
 
 
187
 
 
188
void initroutines(void)
 
189
{
 
190
  Py_InitModule3("routines", module_methods, module_doc);
 
191
  fffpy_import_array();
 
192
  import_array();   /* required NumPy initialization */
 
193
  return; 
 
194
}