← Back to branch summary

~nipy-developers/nipy/fff2

« back to all changes in this revision

Viewing changes to core/fff_gen_stats.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:
core/fff_gen_stats.c
1
 
#include "fff_routines.h"
 
1
#include "fff_gen_stats.h"
 
2
#include "fff_lapack.h"
2
3
 
3
4
#include <string.h>
4
5
#include <stdlib.h>
5
6
#include <math.h>
6
7
#include <errno.h>
7
8
 
8
 
/*
9
 
#include <gsl/gsl_sort_vector.h>
10
 
#include <gsl/gsl_statistics.h>
11
 
#include <gsl/gsl_rng.h>
12
 
#include <gsl/gsl_linalg.h>
13
 
#include <gsl/gsl_blas.h>
14
 
*/
15
 
 
 
9
#if 0 
16
10
 
17
11
extern void fff_rng_draw_noreplace(size_t* x, unsigned int n, unsigned int N)
18
12
{ 
38
32
  return; 
39
33
}
40
34
 
 
35
#endif 
 
36
 
 
37
 
41
38
/* 
42
39
   Squared mahalanobis distance: d2 = x' S^-1 x 
 
40
   Beware: x is not const
43
41
*/ 
44
 
extern double fff_mahalanobis( const fff_vector* x, gsl_matrix* S, fff_vector* work  )
 
42
extern double fff_mahalanobis( fff_vector* x, fff_matrix* S, fff_matrix* Saux  )
45
43
{
46
44
  double d2; 
47
 
 
48
 
  /* Copy x into work */ 
49
 
  fff_vector_memcpy(work, x); 
 
45
  double m = 0.0; 
50
46
 
51
47
  /* Cholesky decomposition: S = L L^t, L lower triangular */
52
 
  gsl_linalg_cholesky_decomp(S); 
 
48
  fff_lapack_dpotrf(CblasLower, S, Saux); 
53
49
 
54
50
  /* Compute S^-1 x */           
55
 
  gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, S, work); /* L^-1 x */ 
56
 
  gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, S, work); /* L^-t (L^-1 x) */ 
 
51
  fff_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, S, x); /* L^-1 x */ 
57
52
  
58
53
  /* Compute x' S^-1 x */ 
59
 
  gsl_blas_ddot (x, work, &d2); 
 
54
  d2 = (double) ( fff_vector_ssd(x, &m, 1) / (long double)x->size ); 
60
55
 
61
 
  return(d2); }
 
56
  return d2; 
 
57
}
62
58