← Back to branch summary

~nipy-developers/nipy/fff2

« back to all changes in this revision

Viewing changes to python/wrapper/fffpy.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/wrapper/fffpy.c
1
 
#include "fffPy.h"
 
1
#include "fffpy.h"
 
2
 
 
3
 
 
4
extern void _no_warning_( void )
 
5
{
 
6
  import_array(); 
 
7
  return; 
 
8
}
 
9
 
2
10
 
3
11
 
4
12
/* This function must be called before the module can work
5
13
   because PyArray_API is defined static, in order not to share that symbol
6
14
   within the dso. (import_array() asks the pointer value to the python process)
7
15
*/
8
 
void fffPy_import_array(void) { 
 
16
void fffpy_import_array(void) { 
9
17
  import_array(); 
 
18
  return;
10
19
}
11
20
 
12
21
/* Static functions */
13
22
static fff_vector* _fff_vector_fromPtr( const char* data, npy_intp dim, npy_intp stride, int type, int isaligned ); 
14
23
static void _fff_vector_copy_fromPtr( fff_vector* y, const char* data, npy_intp stride, int type ); 
15
24
npy_intp _PyArray_find_main_axis( const PyArrayObject* x, int* ok ); 
16
 
static fff_vector* _fff_vector_fromPyArrayIter( const PyArrayIterObject* it, npy_intp axis ); 
17
 
static void _fff_vector_fromPyArrayIter_update( fff_vector* y, const PyArrayIterObject* it, npy_intp axis );
18
 
static void _fff_vector_toPyArrayIter_update( const fff_vector* y, PyArrayIterObject* it, npy_intp axis );
 
25
 
19
26
 
20
27
/* 
21
28
   Get a fff_vector from an input PyArray. This function acts as a
46
53
/* This routine assumes that y->owner is true (and so y is contiguous) */ 
47
54
static void _fff_vector_copy_fromPtr( fff_vector* y, const char* data, npy_intp stride, int type ) 
48
55
{
49
 
  npy_intp dim = (npy_intp) y->size; 
50
 
  PyArrayObject* x = (PyArrayObject*) PyArray_New( &PyArray_Type, 1, &dim, type, &stride, 
51
 
                                                   (void*)data, 0, NPY_BEHAVED, NULL ); 
52
 
  PyArrayObject* xd = (PyArrayObject*) PyArray_SimpleNewFromData( 1, &dim, NPY_DOUBLE, (void*)y->data );
 
56
  npy_intp dim[1] = {(npy_intp)y->size}; 
 
57
  npy_intp strides[1] = {stride};   
 
58
  PyArrayObject* x = (PyArrayObject*) PyArray_New( &PyArray_Type, 1, dim, type, strides, 
 
59
                                                   (void*)data, 0, NPY_CARRAY, NULL ); 
 
60
  PyArrayObject* xd = (PyArrayObject*) PyArray_SimpleNewFromData( 1, dim, NPY_DOUBLE, (void*)y->data );
53
61
  PyArray_CastTo( xd, x ); 
54
62
  Py_XDECREF( xd );
55
63
  Py_XDECREF( x );
256
264
 
257
265
/** Static routines **/ 
258
266
 
259
 
static fff_vector* _fff_vector_fromPyArrayIter( const PyArrayIterObject* it, npy_intp axis )
260
 
{
261
 
  fff_vector* y; 
262
 
  char* data = PyArray_ITER_DATA( it );
263
 
  PyArrayObject* ao = (PyArrayObject*) it->ao; 
264
 
  npy_intp dim = PyArray_DIM( ao, axis ); 
265
 
  npy_intp stride = PyArray_STRIDE( ao, axis ); 
266
 
  int type = PyArray_TYPE( ao ); 
267
 
  int isaligned = PyArray_ISALIGNED( ao ); 
268
 
 
269
 
  y = _fff_vector_fromPtr( data, dim, stride, type, isaligned ); 
270
 
 
271
 
  return y;
272
 
}
273
 
 
274
 
static void _fff_vector_fromPyArrayIter_update( fff_vector* y, const PyArrayIterObject* it, npy_intp axis ) 
275
 
{
276
 
  /* Case owner (PyArray non-DOUBLE): copy data */ 
277
 
  if ( y->owner ) {
278
 
    PyArrayObject* ao = (PyArrayObject*) it->ao; 
279
 
    _fff_vector_copy_fromPtr( y, PyArray_ITER_DATA(it),
280
 
                              PyArray_STRIDE(ao, axis), PyArray_TYPE(ao) ); 
281
 
    
282
 
  }
283
 
  /* Case non-owner (PyArray DOUBLE): borrow data */ 
284
 
  else 
285
 
    y->data = (double*) PyArray_ITER_DATA( it ); 
286
 
 
287
 
  return; 
288
 
}
289
 
 
290
 
/* y->owner assumed true */ 
291
 
static void _fff_vector_toPyArrayIter_update( const fff_vector* y, PyArrayIterObject* it, npy_intp axis ) 
292
 
{
293
 
  npy_intp dims[1], strides[1]; 
294
 
  PyArrayObject* ao = (PyArrayObject*) it->ao; 
295
 
  dims[0] = (npy_intp) y->size;
296
 
  strides[0] = PyArray_STRIDE(ao, axis);
297
 
 
298
 
  /* Encapsulate y into a PyArray */  
299
 
  PyArrayObject* Py_y = (PyArrayObject*) PyArray_SimpleNewFromData( 1, dims, NPY_DOUBLE, (void*)y->data );
300
 
 
301
 
  /* Current view of the PyArray along axis */ 
302
 
  PyArrayObject* Py_view = (PyArrayObject*) PyArray_New( &PyArray_Type, 1, dims, PyArray_TYPE(ao), strides, 
303
 
                                          (void*)PyArray_ITER_DATA(it), 0, NPY_BEHAVED, NULL ); 
304
 
  
305
 
  /* Copy (with cast) Py_y into Py_view */ 
306
 
  PyArray_CastTo( Py_view, Py_y ); 
307
 
 
308
 
  Py_XDECREF( Py_y );
309
 
  Py_XDECREF( Py_view );
310
 
 
311
 
  return; 
312
 
}
313
 
 
314
 
 
315
 
 
316
 
 
317
 
fff_func_iterator* fff_func_iterator_new( npy_intp nargout, PyArrayObject** argout, 
318
 
                                          npy_intp nargin, PyArrayObject** argin,
319
 
                                          npy_intp npy_axis )
320
 
{
321
 
  fff_func_iterator* thisone;
322
 
  npy_intp i; 
323
 
  PyArrayObject *Py_ptr; 
324
 
  int axis = (int)npy_axis; /* For compatibility with the new PyArray_IterAllButAxis API */ 
325
 
 
326
 
  /* Create this one */ 
327
 
  thisone = (fff_func_iterator*)malloc( sizeof(fff_func_iterator) );   
328
 
  if ( thisone == NULL ) 
329
 
    return NULL; 
330
 
  
331
 
  /* Number of args and axis */ 
332
 
  thisone->nargin = nargin; 
333
 
  thisone->nargout = nargout; 
334
 
  thisone->axis = axis; 
335
 
 
336
 
  /* Number of loop items */ 
337
 
  /* TODO: add consistency check -- all arrays should have same dim except in the axis direction */
338
 
  Py_ptr = argin[0]; 
339
 
  thisone->nitems = PyArray_SIZE(Py_ptr) / PyArray_DIM(Py_ptr, axis);
340
 
  
341
 
  /* Allocate local objects (array iterators and local gsl vectors)*/
342
 
  /** Output **/ 
343
 
  thisone->PyIter_argout = (PyArrayIterObject**)malloc( nargout*sizeof(PyArrayIterObject*) ); 
344
 
  thisone->fff_argout = (fff_vector**)malloc( nargout*sizeof(fff_vector*) ); 
345
 
  for( i=0; i<nargout; i++ ) {
346
 
    thisone->PyIter_argout[i] = (PyArrayIterObject*)PyArray_IterAllButAxis( (PyObject*)argout[i], &axis );
347
 
    thisone->fff_argout[i] = _fff_vector_fromPyArrayIter( thisone->PyIter_argout[i], axis );
348
 
  }
349
 
  /** Input arrays **/ 
350
 
  thisone->PyIter_argin = (PyArrayIterObject**)malloc( nargin*sizeof(PyArrayIterObject*) ); 
351
 
  thisone->fff_argin = (fff_vector**)malloc( nargin*sizeof(fff_vector*) ); 
352
 
  for( i=0; i<nargin; i++ ) {
353
 
    thisone->PyIter_argin[i] = (PyArrayIterObject*)PyArray_IterAllButAxis( (PyObject*)argin[i], &axis );
354
 
    thisone->fff_argin[i] = _fff_vector_fromPyArrayIter( thisone->PyIter_argin[i], axis );
355
 
  }
356
 
  
357
 
  return thisone; 
358
 
}
359
 
 
360
 
 
361
 
 
362
 
void fff_func_iterator_delete( fff_func_iterator* thisone )
363
 
{
364
 
  npy_intp i; 
365
 
 
366
 
  for( i=0; i<thisone->nargout; i++ ) {
367
 
    Py_XDECREF( thisone->PyIter_argout[i] ); 
368
 
    fff_vector_delete( thisone->fff_argout[i] ); 
369
 
  }  
370
 
  for( i=0; i<thisone->nargin; i++ ) {
371
 
    Py_XDECREF( thisone->PyIter_argin[i] ); 
372
 
    fff_vector_delete( thisone->fff_argin[i] ); 
373
 
  }  
374
 
  free( thisone->PyIter_argout ); 
375
 
  free( thisone->fff_argout ); 
376
 
  free( thisone->PyIter_argin ); 
377
 
  free( thisone->fff_argin ); 
378
 
  free( thisone ); 
379
 
  
380
 
  return; 
381
 
} 
382
 
 
383
 
 
384
 
 
385
 
 
386
 
void fff_func_iterator_eval( fff_func_iterator* thisone, 
387
 
                             void (*func)( fff_vector**, fff_vector**, void* ), 
388
 
                             void* par )
389
 
{
390
 
  npy_intp i; 
391
 
  npy_intp nargout = thisone->nargout, nargin = thisone->nargin, axis = thisone->axis; 
392
 
  PyArrayIterObject** PyIter_argout = thisone->PyIter_argout; 
393
 
  fff_vector** fff_argout = thisone->fff_argout; 
394
 
  PyArrayIterObject** PyIter_argin = thisone->PyIter_argin; 
395
 
  fff_vector** fff_argin = thisone->fff_argin; 
396
 
  PyArrayIterObject *PyIter_ptr; 
397
 
  fff_vector* fff_ptr;
398
 
 
399
 
  /* Loop over items */ 
400
 
  while( PyIter_argin[0]->index < thisone->nitems ) {
401
 
 
402
 
    /* Do the job for current item */
403
 
    func( fff_argout, fff_argin, par ); 
404
 
 
405
 
    /* Prepare next item */ 
406
 
    for( i=0; i<nargout; i++ ) {
407
 
      PyIter_ptr = PyIter_argout[i];
408
 
      fff_ptr = fff_argout[i]; 
409
 
      
410
 
      /* If the gsl vector is not a view (owner), then push back its content */
411
 
      if ( fff_ptr->owner )     
412
 
        _fff_vector_toPyArrayIter_update( fff_ptr, PyIter_ptr, axis ); 
413
 
      
414
 
      PyArray_ITER_NEXT( PyIter_ptr );
415
 
      _fff_vector_fromPyArrayIter_update( fff_ptr, PyIter_ptr, axis ); 
416
 
    }
417
 
    for( i=0; i<nargin; i++ ) {
418
 
      PyIter_ptr = PyIter_argin[i];
419
 
      PyArray_ITER_NEXT( PyIter_ptr );
420
 
      _fff_vector_fromPyArrayIter_update( fff_argin[i], PyIter_ptr, axis ); 
421
 
    }
422
 
    
423
 
  } /* End loop over items */ 
424
 
  
425
 
  return;
426
 
}
427
 
 
428
 
 
429
 
 
430
 
void fff_func_iterator_reset( fff_func_iterator* thisone )
431
 
{
432
 
  npy_intp i; 
433
 
  npy_intp nargout = thisone->nargout, nargin = thisone->nargin, axis = thisone->axis; 
434
 
  PyArrayIterObject** PyIter_argout = thisone->PyIter_argout; 
435
 
  fff_vector** fff_argout = thisone->fff_argout; 
436
 
  PyArrayIterObject** PyIter_argin = thisone->PyIter_argin; 
437
 
  fff_vector** fff_argin = thisone->fff_argin; 
438
 
  PyArrayIterObject *PyIter_ptr; 
439
 
  fff_vector *fff_ptr; 
440
 
 
441
 
  for( i=0; i<nargout; i++ ) {
442
 
    PyIter_ptr = PyIter_argout[i];
443
 
    fff_ptr = fff_argout[i]; 
444
 
    PyArray_ITER_RESET( PyIter_ptr );
445
 
    _fff_vector_fromPyArrayIter_update( fff_ptr, PyIter_ptr, axis ); 
446
 
  }
447
 
  for( i=0; i<nargin; i++ ) {
448
 
    PyIter_ptr = PyIter_argin[i];
449
 
    fff_ptr = fff_argin[i]; 
450
 
    PyArray_ITER_RESET( PyIter_ptr );
451
 
    _fff_vector_fromPyArrayIter_update( fff_ptr, PyIter_ptr, axis ); 
452
 
  }
453
 
  
454
 
  return; 
455
 
}
 
267
 
456
268
 
457
269
 
458
270
 
633
445
 
634
446
  return x;
635
447
}
 
448
 
 
449
 
 
450
 
 
451
 
 
452
/* 
 
453
   Create a PyArrayMultiArrayIter so as to iterate simultaneously on
 
454
   several arrays except in one common axis.  This is a
 
455
   PyArrayMultiArrayIter constructor that I would like to see in the
 
456
   numpy C API, but for now it is "hosted" in fffPy.h
 
457
 
 
458
   TODO: check dimensions. 
 
459
*/ 
 
460
#include <stdarg.h>
 
461
 
 
462
PyObject* PyArray_MultiIterAllButAxis( int n, int axis, ... )
 
463
{
 
464
  va_list va;
 
465
  PyArrayMultiIterObject *multi;
 
466
  PyObject *current;
 
467
  PyObject *arr;
 
468
  
 
469
  int i, err=0;
 
470
  
 
471
  if (n < 2 || n > NPY_MAXARGS) {
 
472
    PyErr_Format(PyExc_ValueError,
 
473
                 "Need between 2 and (%d) "                             \
 
474
                 "array objects (inclusive).", NPY_MAXARGS);
 
475
    return NULL;
 
476
  }
 
477
  
 
478
  /* fprintf(stderr, "multi new...");*/
 
479
  
 
480
  multi = (PyArrayMultiIterObject *) malloc(sizeof(PyArrayMultiIterObject));
 
481
  if (multi == NULL) return PyErr_NoMemory();
 
482
  PyObject_Init((PyObject *)multi, &PyArrayMultiIter_Type);
 
483
  
 
484
  for (i=0; i<n; i++) multi->iters[i] = NULL;
 
485
  multi->numiter = n;
 
486
  multi->index = 0;
 
487
  
 
488
  va_start(va, axis);
 
489
  for (i=0; i<n; i++) {
 
490
    current = va_arg(va, PyObject *);
 
491
    arr = PyArray_FROM_O(current);
 
492
    if (arr==NULL) {
 
493
      err=1; break;
 
494
    }
 
495
    else {
 
496
      multi->iters[i] = (PyArrayIterObject *)PyArray_IterAllButAxis(arr, &axis);
 
497
      Py_DECREF(arr);
 
498
    }
 
499
  }
 
500
  
 
501
  va_end(va);
 
502
  
 
503
  /*if (!err && PyArray_Broadcast(multi) < 0) err=1;*/ 
 
504
  
 
505
  if (err) {
 
506
    Py_DECREF(multi);
 
507
    return NULL;
 
508
  }
 
509
  
 
510
  PyArray_MultiIter_RESET(multi);
 
511
  
 
512
  return (PyObject *)multi;
 
513
}
 
514
 
 
515
 
 
516
/* Create an fff_vector from a PyArrayIter object */ 
 
517
fff_vector* fff_vector_fromPyArrayIter( const PyArrayIterObject* it, npy_intp axis )
 
518
{
 
519
  fff_vector* y; 
 
520
  char* data = PyArray_ITER_DATA( it );
 
521
  PyArrayObject* ao = (PyArrayObject*) it->ao; 
 
522
  npy_intp dim = PyArray_DIM( ao, axis ); 
 
523
  npy_intp stride = PyArray_STRIDE( ao, axis ); 
 
524
  int type = PyArray_TYPE( ao ); 
 
525
  int isaligned = PyArray_ISALIGNED( ao ); 
 
526
 
 
527
  y = _fff_vector_fromPtr( data, dim, stride, type, isaligned ); 
 
528
 
 
529
  return y;
 
530
}
 
531
 
 
532
/* Fetch vector data in an iterator */ 
 
533
void fff_vector_get_fromPyArrayIter( fff_vector* y, const PyArrayIterObject* it, npy_intp axis ) 
 
534
{
 
535
  /* Case owner (PyArray non-DOUBLE): copy data */ 
 
536
  if ( y->owner ) {
 
537
    PyArrayObject* ao = (PyArrayObject*) it->ao; 
 
538
    _fff_vector_copy_fromPtr( y, PyArray_ITER_DATA(it),
 
539
                              PyArray_STRIDE(ao, axis), PyArray_TYPE(ao) ); 
 
540
    
 
541
  }
 
542
  /* Case non-owner (PyArray DOUBLE): borrow data */ 
 
543
  else 
 
544
    y->data = (double*) PyArray_ITER_DATA( it ); 
 
545
 
 
546
  return; 
 
547
}
 
548
 
 
549
/* y assumed owner (otherwise, there is no point to call this function) */ 
 
550
void fff_vector_const_toPyArrayIter( const fff_vector* y, PyArrayIterObject* it, npy_intp axis ) 
 
551
{
 
552
  if ( ! y->owner )
 
553
    return; 
 
554
 
 
555
  npy_intp dims[1], strides[1]; 
 
556
  PyArrayObject* ao = (PyArrayObject*) it->ao; 
 
557
  dims[0] = (npy_intp) y->size;
 
558
  strides[0] = PyArray_STRIDE(ao, axis);
 
559
 
 
560
  /* Encapsulate y into a PyArray */  
 
561
  PyArrayObject* Py_y = (PyArrayObject*) PyArray_SimpleNewFromData( 1, dims, NPY_DOUBLE, (void*)y->data );
 
562
 
 
563
  /* Current view of the PyArray along axis */ 
 
564
  PyArrayObject* Py_view = (PyArrayObject*) PyArray_New( &PyArray_Type, 1, dims, PyArray_TYPE(ao), strides, 
 
565
                                                         (void*)PyArray_ITER_DATA(it), 0, NPY_BEHAVED, NULL ); 
 
566
  
 
567
  /* Copy (with cast) Py_y into Py_view */ 
 
568
  PyArray_CastTo( Py_view, Py_y ); 
 
569
 
 
570
  Py_XDECREF( Py_y );
 
571
  Py_XDECREF( Py_view );
 
572
 
 
573
  return; 
 
574
}
 
575
 
 
576
 
 
577