MNE-CPP  0.1.9
A Framework for Electrophysiology
fiff_sparse_matrix.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "fiff_sparse_matrix.h"
42 #include <fiff/fiff_file.h>
43 #include <fiff/fiff_types.h>
44 
45 #define MALLOC_18(x,t) (t *)malloc((x)*sizeof(t))
46 
47 #define REALLOC_18(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
48 
49 #define FREE_18(x) if ((char *)(x) != NULL) free((char *)(x))
50 
51 //=============================================================================================================
52 // USED NAMESPACES
53 //=============================================================================================================
54 
55 using namespace Eigen;
56 using namespace FIFFLIB;
57 using namespace FIFFLIB;
58 
59 //============================= fiff_type_spec.h =============================
60 
61 /*
62  * These return information about a fiff type.
63  */
64 
65 fiff_int_t fiff_type_base(fiff_int_t type)
66 {
67  return type & FIFFTS_BASE_MASK;
68 }
69 
70 fiff_int_t fiff_type_fundamental(fiff_int_t type)
71 {
72  return type & FIFFTS_FS_MASK;
73 }
74 
75 fiff_int_t fiff_type_matrix_coding(fiff_int_t type)
76 {
77  return type & FIFFTS_MC_MASK;
78 }
79 
80 //============================= fiff_matrix.c =============================
81 
82 int *fiff_get_matrix_dims(FiffTag::SPtr& tag)
83 /*
84  * Interpret dimensions from matrix data (dense and sparse)
85  */
86 {
87  int ndim;
88  int *dims;
89  int *res,k;
90  unsigned int tsize = tag->size();
91  /*
92  * Initial checks
93  */
94  if (tag->data() == NULL) {
95  qCritical("fiff_get_matrix_dims: no data available!");
96  return NULL;
97  }
98  if (fiff_type_fundamental(tag->type) != FIFFTS_FS_MATRIX) {
99  qCritical("fiff_get_matrix_dims: tag does not contain a matrix!");
100  return NULL;
101  }
102  if (tsize < sizeof(fiff_int_t)) {
103  qCritical("fiff_get_matrix_dims: too small matrix data!");
104  return NULL;
105  }
106  /*
107  * Get the number of dimensions and check
108  */
109  ndim = *((fiff_int_t *)((fiff_byte_t *)(tag->data())+tag->size()-sizeof(fiff_int_t)));
110  if (ndim <= 0 || ndim > FIFFC_MATRIX_MAX_DIM) {
111  qCritical("fiff_get_matrix_dims: unreasonable # of dimensions!");
112  return NULL;
113  }
114  if (fiff_type_matrix_coding(tag->type) == FIFFTS_MC_DENSE) {
115  if (tsize < (ndim+1)*sizeof(fiff_int_t)) {
116  qCritical("fiff_get_matrix_dims: too small matrix data!");
117  return NULL;
118  }
119  res = MALLOC_18(ndim+1,int);
120  res[0] = ndim;
121  dims = ((fiff_int_t *)((fiff_byte_t *)(tag->data())+tag->size())) - ndim - 1;
122  for (k = 0; k < ndim; k++)
123  res[k+1] = dims[k];
124  }
125  else if (fiff_type_matrix_coding(tag->type) == FIFFTS_MC_CCS ||
126  fiff_type_matrix_coding(tag->type) == FIFFTS_MC_RCS) {
127  if (tsize < (ndim+2)*sizeof(fiff_int_t)) {
128  qCritical("fiff_get_matrix_sparse_dims: too small matrix data!");
129  return NULL; }
130 
131  res = MALLOC_18(ndim+2,int);
132  res[0] = ndim;
133  dims = ((fiff_int_t *)((fiff_byte_t *)(tag->data())+tag->size())) - ndim - 1;
134  for (k = 0; k < ndim; k++)
135  res[k+1] = dims[k];
136  res[ndim+1] = dims[-1];
137  }
138  else {
139  qCritical("fiff_get_matrix_dims: unknown matrix coding.");
140  return NULL;
141  }
142  return res;
143 }
144 
145 //=============================================================================================================
146 // DEFINE MEMBER METHODS
147 //=============================================================================================================
148 
149 FiffSparseMatrix::FiffSparseMatrix()
150 : coding(0)
151 , m(0)
152 , n(0)
153 , nz(0)
154 , data(NULL)
155 , inds(NULL)
156 , ptrs(NULL)
157 {
158 }
159 
160 //=============================================================================================================
161 
163 {
164  int size;
165 
166  this->coding = mat.coding;
167  this->m = mat.m;
168  this->n = mat.n;
169  this->nz = mat.nz;
170 
171  if (mat.coding == FIFFTS_MC_CCS) {
172  size = mat.nz*(sizeof(FIFFLIB::fiff_float_t) + sizeof(FIFFLIB::fiff_int_t)) +
173  (mat.n+1)*(sizeof(FIFFLIB::fiff_int_t));
174  }
175  if (mat.coding == FIFFTS_MC_RCS) {
176  size = mat.nz*(sizeof(FIFFLIB::fiff_float_t) + sizeof(FIFFLIB::fiff_int_t)) +
177  (mat.m+1)*(sizeof(FIFFLIB::fiff_int_t));
178  }
179  else {
180  printf("Illegal sparse matrix storage type: %d",mat.coding);
181  return;
182  }
183  this->data = (float *)malloc(size);
184  this->inds = (int *)(this->data+this->nz);
185  this->ptrs = this->inds+this->nz;
186  memcpy(data,mat.data,size);
187 }
188 
189 //=============================================================================================================
190 
192 {
193  if(data)
194  FREE_18(data);
195 }
196 
197 //=============================================================================================================
198 
199 fiff_int_t *FiffSparseMatrix::fiff_get_matrix_sparse_dims(FiffTag::SPtr &tag)
200 {
201  return fiff_get_matrix_dims(tag);
202 }
203 
204 //=============================================================================================================
205 
206 FiffSparseMatrix *FiffSparseMatrix::fiff_get_float_sparse_matrix(FiffTag::SPtr &tag)
207 {
208  int *dims;
209  FIFFLIB::FiffSparseMatrix* res = NULL;
210  int m,n,nz;
211  int coding,correct_size;
212 
213  if ( fiff_type_fundamental(tag->type) != FIFFT_MATRIX ||
214  fiff_type_base(tag->type) != FIFFT_FLOAT ||
215  (fiff_type_matrix_coding(tag->type) != FIFFTS_MC_CCS &&
216  fiff_type_matrix_coding(tag->type) != FIFFTS_MC_RCS) ) {
217  printf("fiff_get_float_ccs_matrix: wrong data type!");
218  return NULL;
219  }
220 
221  if ((dims = fiff_get_matrix_sparse_dims(tag)) == NULL)
222  return NULL;
223 
224  if (dims[0] != 2) {
225  printf("fiff_get_float_sparse_matrix: wrong # of dimensions!");
226  return NULL;
227  }
228 
229  m = dims[1];
230  n = dims[2];
231  nz = dims[3];
232 
233  coding = fiff_type_matrix_coding(tag->type);
234  if (coding == FIFFTS_MC_CCS)
235  correct_size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
236  (n+1+dims[0]+2)*(sizeof(fiff_int_t));
237  else if (coding == FIFFTS_MC_RCS)
238  correct_size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
239  (m+1+dims[0]+2)*(sizeof(fiff_int_t));
240  else {
241  printf("fiff_get_float_sparse_matrix: Incomprehensible sparse matrix coding");
242  return NULL;
243  }
244  if (tag->size() != correct_size) {
245  printf("fiff_get_float_sparse_matrix: wrong data size!");
246  FREE_18(dims);
247  return NULL;
248  }
249  /*
250  * Set up structure
251  */
252  res = new FIFFLIB::FiffSparseMatrix;
253  res->m = m;
254  res->n = n;
255  res->nz = nz;
256  res->data = MALLOC_18(correct_size,float);
257  memcpy (res->data,(float*)tag->data(),correct_size);
258  res->coding = coding;
259  res->inds = (int *)(res->data + res->nz);
260  res->ptrs = res->inds + res->nz;
261 
262  FREE_18(dims);
263 
264  return res;
265 }
266 
267 //=============================================================================================================
268 
269 FiffSparseMatrix *FiffSparseMatrix::create_sparse_rcs(int nrow, int ncol, int *nnz, int **colindex, float **vals) /* The nonzero elements on each row
270  * If null, the matrix will be all zeroes */
271 {
272  FIFFLIB::FiffSparseMatrix* sparse = NULL;
273  int j,k,nz,ptr,size,ind;
274  int stor_type = FIFFTS_MC_RCS;
275 
276  for (j = 0, nz = 0; j < nrow; j++)
277  nz = nz + nnz[j];
278 
279  if (nz <= 0) {
280  printf("No nonzero elements specified.");
281  return NULL;
282  }
283  if (stor_type == FIFFTS_MC_RCS) {
284  size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
285  (nrow+1)*(sizeof(fiff_int_t));
286  }
287  else {
288  printf("Illegal sparse matrix storage type: %d",stor_type);
289  return NULL;
290  }
291  sparse = new FIFFLIB::FiffSparseMatrix;
292  sparse->coding = stor_type;
293  sparse->m = nrow;
294  sparse->n = ncol;
295  sparse->nz = nz;
296  sparse->data = (float *)malloc(size);
297  sparse->inds = (int *)(sparse->data+nz);
298  sparse->ptrs = sparse->inds+nz;
299 
300  for (j = 0, nz = 0; j < nrow; j++) {
301  ptr = -1;
302  for (k = 0; k < nnz[j]; k++) {
303  if (ptr < 0)
304  ptr = nz;
305  ind = sparse->inds[nz] = colindex[j][k];
306  if (ind < 0 || ind >= ncol) {
307  printf("Column index out of range in mne_create_sparse_rcs");
308  goto bad;
309  }
310  if (vals)
311  sparse->data[nz] = vals[j][k];
312  else
313  sparse->data[nz] = 0.0;
314  nz++;
315  }
316  sparse->ptrs[j] = ptr;
317  }
318  sparse->ptrs[nrow] = nz;
319  for (j = nrow-1; j >= 0; j--) /* Take care of the empty rows */
320  if (sparse->ptrs[j] < 0)
321  sparse->ptrs[j] = sparse->ptrs[j+1];
322  return sparse;
323 
324 bad : {
325  if(sparse)
326  delete sparse;
327  return NULL;
328  }
329 }
330 
331 //=============================================================================================================
332 
333 FiffSparseMatrix *FiffSparseMatrix::mne_add_upper_triangle_rcs()
334 /*
335  * Fill in upper triangle with the lower triangle values
336  */
337 {
338  int *nnz = NULL;
339  int **colindex = NULL;
340  float **vals = NULL;
341  FIFFLIB::FiffSparseMatrix* res = NULL;
342  int i,j,k,row;
343  int *nadd = NULL;
344 
345  if (this->coding != FIFFTS_MC_RCS) {
346  printf("The input matrix to mne_add_upper_triangle_rcs must be in RCS format");
347  goto out;
348  }
349  if (this->m != this->n) {
350  printf("The input matrix to mne_add_upper_triangle_rcs must be square");
351  goto out;
352  }
353  nnz = MALLOC_18(this->m,int);
354  colindex = MALLOC_18(this->m,int *);
355  vals = MALLOC_18(this->m,float *);
356  for (i = 0; i < this->m; i++) {
357  nnz[i] = this->ptrs[i+1] - this->ptrs[i];
358  if (nnz[i] > 0) {
359  colindex[i] = MALLOC_18(nnz[i],int);
360  vals[i] = MALLOC_18(nnz[i],float);
361  for (j = this->ptrs[i], k = 0; j < this->ptrs[i+1]; j++, k++) {
362  vals[i][k] = this->data[j];
363  colindex[i][k] = this->inds[j];
364  }
365  }
366  else {
367  colindex[i] = NULL;
368  vals[i] = NULL;
369  }
370  }
371  /*
372  * Add the elements
373  */
374  nadd = MALLOC_18(this->m,int);
375  for (i = 0; i < this->m; i++)
376  nadd[i] = 0;
377  for (i = 0; i < this->m; i++)
378  for (j = this->ptrs[i]; j < this->ptrs[i+1]; j++)
379  nadd[this->inds[j]]++;
380  for (i = 0; i < this->m; i++) {
381  colindex[i] = REALLOC_18(colindex[i],nnz[i]+nadd[i],int);
382  vals[i] = REALLOC_18(vals[i],nnz[i]+nadd[i],float);
383  }
384  for (i = 0; i < this->m; i++)
385  for (j = this->ptrs[i]; j < this->ptrs[i+1]; j++) {
386  row = this->inds[j];
387  colindex[row][nnz[row]] = i;
388  vals[row][nnz[row]] = this->data[j];
389  nnz[row]++;
390  }
391  res = create_sparse_rcs(this->m,this->n,nnz,colindex,vals);
392 
393 out : {
394  for (i = 0; i < this->m; i++) {
395  FREE_18(colindex[i]);
396  FREE_18(vals[i]);
397  }
398  FREE_18(nnz);
399  FREE_18(vals);
400  FREE_18(colindex);
401  FREE_18(nadd);
402  return res;
403  }
404 }
Old fiff_type declarations - replace them.
Header file describing the numerical values used in fif files.
FIFFLIB::fiff_int_t coding
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
FiffSparseMatrix class declaration.
FIFFLIB::fiff_int_t * ptrs
FIFFLIB::fiff_int_t * inds
FIFFLIB::fiff_float_t * data
Data associated with MNE computations for each mneMeasDataSet.