← Back to branch summary

~nipy-developers/nipy/fff2

« back to all changes in this revision

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