← Back to branch summary

~nipy-developers/nipy/fff2

« back to all changes in this revision

Viewing changes to matlab/fffMat_onesample_stat.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:
matlab/fffMat_onesample_stat.c
 
1
#include "fffMat.h"
 
2
 
 
3
#include <fff_onesample_stat.h>
 
4
#include <gsl/gsl_rng.h>
 
5
 
 
6
#include <string.h>
 
7
#include <time.h>
 
8
 
 
9
#define DEFAULT_FLAG FFF_ONESAMPLE_STUDENT
 
10
#define DEFAULT_BASE 0.0
 
11
#define DEFAULT_NSIMU 1
 
12
#define DEFAULT_MAGIC 0
 
13
#define DEFAULT_NITER 5
 
14
 
 
15
#define MFX_OFFSET 10 
 
16
 
 
17
/* 
 
18
   
 
19
T = fffMat_onesample_stat( E, V, flagStat, base, nsimu, magic, niter ); 
 
20
 
 
21
magic < 0    ==> randomize permutations
 
22
magic >= 0   ==> produce deterministic permutations in the range [magic, magic+nsimu[
 
23
 
 
24
*/ 
 
25
 
 
26
static void _set_simulation_parameters( int* flag_permute, int* flag_use_rng, 
 
27
                                        int* simu0, int* nsimu, 
 
28
                                        size_t* nsimu_max, int n, 
 
29
                                        long int magic )
 
30
{
 
31
  if ( *nsimu < 1 )
 
32
    *nsimu = 1;     
 
33
  *flag_permute = ( (*nsimu > 1) || (magic!=0) ); 
 
34
  *flag_use_rng = ( magic < 0 );
 
35
  
 
36
  *simu0 = 0; 
 
37
  if ( magic == 0 )
 
38
    *simu0 = 1; 
 
39
 
 
40
  *nsimu_max = 1 << n; 
 
41
  if ( *flag_permute ) 
 
42
    if ( ( (size_t)(*nsimu) > *nsimu_max ) && 
 
43
         ( !(*flag_use_rng) ) )
 
44
      *nsimu = *nsimu_max;
 
45
  
 
46
  return; 
 
47
}
 
48
 
 
49
void mexFunction( int nlhs, mxArray *plhs[], 
 
50
                  int nrhs, const mxArray*prhs[] )
 
51
{ 
 
52
  /* Variables */ 
 
53
  int flag = DEFAULT_FLAG;
 
54
  double base = DEFAULT_BASE; 
 
55
  int nsimu = DEFAULT_NSIMU;
 
56
  long int magic = DEFAULT_MAGIC; 
 
57
  int niter = DEFAULT_NITER;
 
58
  gsl_vector **E, **V, **T;
 
59
  int vox, nvox, nV, simu, simu0;
 
60
  fff_onesample_stat* Stat; 
 
61
  fff_onesample_stat_mfx* StatMFX; 
 
62
  size_t nsimu_max; 
 
63
  int flag_permute, flag_use_rng, flag_mfx = 0;
 
64
  long int magic_aux; 
 
65
  gsl_rng* rng = NULL;
 
66
  double *aux; 
 
67
 
 
68
  /* Check for proper number of arguments. */
 
69
  if ( nrhs < 2 ) {
 
70
    mexErrMsgTxt("Two input arguments required.");
 
71
  } 
 
72
  else if ( nlhs > 50 ) {
 
73
    mexErrMsgTxt("Too many output arguments");
 
74
  }
 
75
 
 
76
  /* Get list of effect vectors */
 
77
  E = gsl_vector_list_from_mxArray_new( prhs[0], &nvox );
 
78
  
 
79
  /* Get list of associated variances */
 
80
  V = gsl_vector_list_from_mxArray_new( prhs[1], &nV );
 
81
  
 
82
  /* Test consistency between E and V */
 
83
  if ( (nV != 0) && (nV != nvox) ) {
 
84
    gsl_vector_list_from_mxArray_delete( E, nvox );
 
85
    gsl_vector_list_from_mxArray_delete( V, nV );
 
86
    mexErrMsgTxt("Incompatible effect and variance arrays.");
 
87
  }
 
88
 
 
89
  /* Statistic choice */
 
90
  if ( nrhs > 2 ) {
 
91
    aux = mxGetPr( prhs[2] );
 
92
    flag = (int)(*aux);
 
93
    if ( flag >= MFX_OFFSET ) 
 
94
      flag_mfx = 1; 
 
95
  }
 
96
  
 
97
  /* Baseline */
 
98
  if ( nrhs > 3 ) {
 
99
    aux = mxGetPr( prhs[3] );
 
100
    base = (double)(*aux);
 
101
  }
 
102
 
 
103
  /* Number of simulations */
 
104
  if ( nrhs > 4 ) {
 
105
    aux = mxGetPr( prhs[4] );
 
106
    nsimu = (int)(*aux);
 
107
  }
 
108
 
 
109
  /* Magic number */ 
 
110
  if ( nrhs > 5 ) {
 
111
    aux = mxGetPr( prhs[5] );
 
112
    magic = (long int)(*aux);
 
113
  }
 
114
 
 
115
  /* Number of EM iterations (MFX stats) */ 
 
116
  if ( nrhs > 6 ) {
 
117
    aux = mxGetPr( prhs[6] );
 
118
    niter = (int)(*aux);
 
119
  }
 
120
 
 
121
  
 
122
  /* Simulation parameters */ 
 
123
  _set_simulation_parameters( &flag_permute, &flag_use_rng, 
 
124
                              &simu0, &nsimu, &nsimu_max, 
 
125
                              (int)E[0]->size, magic ); 
 
126
 
 
127
  /* Output */ 
 
128
  plhs[0] = mxCreateDoubleMatrix( nsimu, nvox, mxREAL );
 
129
  
 
130
  /* Local structures */
 
131
  T = gsl_vector_list_from_mxArray_new( plhs[0], &nvox );
 
132
  if ( ! flag_mfx )
 
133
    Stat = fff_onesample_stat_new( (size_t)E[0]->size, (fff_onesample_stat_flag)flag, base );
 
134
  else {
 
135
    StatMFX = fff_onesample_stat_mfx_new( (size_t)E[0]->size, (fff_onesample_stat_mfx_flag)(flag-MFX_OFFSET), base );
 
136
    fff_onesample_stat_mfx_set_niter( StatMFX, (unsigned int)niter );
 
137
  }
 
138
  if ( flag_use_rng ) {
 
139
    gsl_rng_env_setup();
 
140
    rng = gsl_rng_alloc( gsl_rng_default );
 
141
  }
 
142
 
 
143
  /* Compute stat for the original dataset */ 
 
144
  if ( simu0 == 1 ) 
 
145
    if ( ! flag_mfx ) 
 
146
      for ( vox=0; vox<nvox; vox++ ) {
 
147
        fff_onesample_stat_eval( Stat, E[vox] );
 
148
        T[vox]->data[0] = Stat->stat; 
 
149
      }
 
150
    else 
 
151
      for ( vox=0; vox<nvox; vox++ ) {
 
152
        fff_onesample_stat_mfx_eval( StatMFX, E[vox], V[vox] );
 
153
        T[vox]->data[0] = StatMFX->stat; 
 
154
      }
 
155
 
 
156
  /* Simulations */ 
 
157
  for ( simu=simu0; simu<nsimu; simu++ ) {
 
158
 
 
159
    /* Determine permutation (magic number) */ 
 
160
    if ( flag_use_rng ) 
 
161
      magic_aux = gsl_rng_uniform_int( rng, nsimu_max );
 
162
    else 
 
163
      magic_aux = (size_t)(simu + magic); 
 
164
    
 
165
    /* Loop over voxels */
 
166
    if ( ! flag_mfx ) {
 
167
      Stat->magic = magic_aux; 
 
168
      for ( vox=0; vox<nvox; vox++ ) {
 
169
        fff_onesample_stat_permute_and_eval( Stat, E[vox] );   
 
170
        T[vox]->data[simu] = Stat->stat;
 
171
      }
 
172
    }
 
173
    else {
 
174
      StatMFX->magic = magic_aux; 
 
175
      for ( vox=0; vox<nvox; vox++ ) {
 
176
        fff_onesample_stat_mfx_permute_and_eval( StatMFX, E[vox], V[vox] );   
 
177
        T[vox]->data[simu] = StatMFX->stat;
 
178
      }
 
179
    }
 
180
  
 
181
  } /* End loop over simulations */ 
 
182
 
 
183
  /* Free memory */ 
 
184
  gsl_vector_list_from_mxArray_delete( E, nvox );
 
185
  gsl_vector_list_from_mxArray_delete( V, nV );
 
186
  gsl_vector_list_from_mxArray_delete( T, nvox );
 
187
  if ( ! flag_mfx ) 
 
188
    fff_onesample_stat_delete( Stat ); 
 
189
  else 
 
190
    fff_onesample_stat_mfx_delete( StatMFX ); 
 
191
  if ( flag_use_rng ) 
 
192
    gsl_rng_free( rng ); 
 
193
 
 
194
  return; 
 
195
}
 
196