src/averaging.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2008 by Matthias Ihrke   *
 *   mihrke@uni-goettingen.de   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "averaging.h"
#include "clustering.h"
#include "eeg.h"
#include "linalg.h"

Array* average_example( const Array *data, uint idx[2], double weights[2], OptArgList *optargs ){
  int i,j;
  if( data->ndim!=3 || data->dtype!=DOUBLE ){
     errprintf("data must be 3D and DOUBLE, have %i, %i\n", data->ndim, data->dtype );
     return NULL;
  } 
  int N,C,n;
  N = data->size[0];
  C = data->size[1];
  n = data->size[2];
  Array *avg = array_new2( DOUBLE, 2, C, n );
  for( i=0; i<C; i++ ){
     for( j=0; j<n; j++ ){
        mat_IDX( avg, i, j )=
          (weights[0]*array_INDEX3( data, double, idx[0], i, j ))+
          (weights[1]*array_INDEX3( data, double, idx[1], i, j ));
     }
  }                                 

  return avg;
}
Array* average_warpmarkers( const Array *data, uint idx[2], double weights[2], OptArgList *optargs ){
}

Array* hierarchical_average( const Array *data, const Array *distmat, 
                                      SignalAverageFunction avgfct, OptArgList *optargs ){
  Dendrogram *T, *Tsub; 
  ProgressBarFunction progress=NULL;  
  LinkageFunction linkage=dgram_dist_completelinkage;
  void *ptr;

  bool ismatrix;
  matrix_CHECKSQR( ismatrix, distmat );
  if( !ismatrix ) return NULL;
  if( data->ndim>3 || data->ndim<2 || data->dtype!=DOUBLE ){
     errprintf( "Data needs to be 2D or 3D, got %iD and of type DOUBLE\n", data->ndim );
     return NULL;
  }

  optarg_PARSE_PTR( optargs, "linkage", linkage, LinkageFunction, ptr );
  optarg_PARSE_PTR( optargs, "progress", progress, ProgressBarFunction, ptr );

  int N,C,n;
  N = data->size[0];
  C = (data->ndim==2)?1:(data->size[1]);
  n = (data->ndim==2)?(data->size[1]):(data->size[2]);
  dprintf("Data is (%i x %i x %i)\n", N, C, n );

  /* wrapper for 2D-data -> work on d instead of data */
  Array *d=array_copy( data, TRUE );
  if( d->ndim==2 ){
     d->ndim=3;
     free( d->size );
     d->size=(uint*)malloc( 3*sizeof(uint) );
     d->size[0]=N; d->size[1]=C; d->size[2]=n;
  }

  /*------------------- computation --------------------------*/
  /* build hierarchical dendrogram */
  T = agglomerative_clustering( distmat, linkage );

  int i;
  double weights[2];
  int *indices = (int*)malloc( N*sizeof(int) );
  for( i=0; i<N; i++ ){
     indices[i] = 1;
  }  

  /* now walk the tree to find pairs of trials to match */
  Array *new, 
     *avg=array_new2( DOUBLE, 2, C, n );
  int trials_left = N;
  uint idx[2]={0,0};
  if( progress ) progress( PROGRESSBAR_INIT, N );

  while( trials_left >= 2 ){
     if( progress ) progress( PROGRESSBAR_CONTINUE_LONG, N-trials_left );
     
     Tsub = dgram_get_deepest( T );
     idx[0] = Tsub->left->val;
     idx[1] = Tsub->right->val;
     
     weights[0] = indices[idx[0]]/(double)(indices[idx[0]]+indices[idx[1]]);
     weights[1] = indices[idx[1]]/(double)(indices[idx[0]]+indices[idx[1]]);
     /* Array* average_example( const AADTWrray *data, uint idx[2], 
               double weights[2], OptArgList *optargs ) */
     new=avgfct( d, idx, weights, optargs );
     memcpy( array_INDEXMEM3( d, idx[0], 0, 0 ),
                new->data, C*n*sizeof( double ) );
     array_free( new );

     indices[idx[0]] += indices[idx[1]];
     indices[idx[1]]=-1;              /* never to be used again */
     
     /* replace node by leaf representing average(idx1, idx2) */
     Tsub->val = idx[0];
     Tsub->left = NULL;
     Tsub->right = NULL;    
     trials_left--;
  }
  /* idx[0] is the final average */
  memcpy( avg->data, array_INDEXMEM3( d, idx[0], 0, 0 ),
             C*n*sizeof( double ) );
  /*------------------- /computation --------------------------*/

  free( indices );
  array_free( d );

  array_dimred( avg );
  return avg;
}



EEG*     eeg_simple_average( const EEG *eeg ){
#ifdef FIXEEG
  EEG *out;
  int c;

  out = eeg_init( eeg->nbchan, 1, eeg->n ); /* the average */
  out->sampling_rate=eeg->sampling_rate;
  if( eeg->times ){
     out->times = (double*) malloc( eeg->n*sizeof(double) );
     memcpy( out->times, eeg->times, eeg->n*sizeof(double) );
  }
  if( eeg->chaninfo ){
     out->chaninfo = (ChannelInfo*) malloc( eeg->nbchan*sizeof(ChannelInfo) );
     memcpy( out->chaninfo, eeg->chaninfo,  eeg->nbchan*sizeof(ChannelInfo) );
  }
  eeg_append_comment( out, "output from eegtrials_simple_average()\n" );
  for( c=0; c<eeg->nbchan; c++ ){
     /* TODO */
  }
  return out;
#endif
}


#ifdef EXPERIMENTAL

EEG*     eeg_average_channels( const EEG *eeg ){
#ifdef FIXEEG

  EEG *out;
  int c, i, j;

  out = eeg_init( 1, eeg->ntrials, eeg->n ); /* the average */
  out->sampling_rate=eeg->sampling_rate;
  if( eeg->times ){
     out->times = (double*) malloc( eeg->n*sizeof(double) );
     memcpy( out->times, eeg->times, eeg->n*sizeof(double) );
  }
  eeg_append_comment( out, "output from eegtrials_average_channels()\n" );

  for( i=0; i<eeg->ntrials; i++ ){
     for( j=0; j<eeg->n; j++ ){
        out->data[0][i][j] = 0;
        for( c=0; c<eeg->nbchan; c++ ){
          out->data[0][i][j] += eeg->data[c][i][j];
        }
        out->data[0][i][j] /= (double)eeg->nbchan;
     }
  }
  return out;
#endif
}

Array* average_unequal_warp( Array **data, uint idx[2], double weights[2], OptArgList *optargs ){
  int i,j;

  int N,C,n;
  C = data[idx[0]]->size[1];
  Array *a=data[idx[0]];
  Array *b=data[idx[1]];

  Array *d = distmatrix_signaldist( vectordist_euclidean, 
                                                a, b, NULL, NULL );

  OptArgList *opts=optarglist( "slope_constraint=int", SLOPE_CONSTRAINT_NONE );
  Array *m1 = matrix_dtw_cumulate( d, TRUE, opts );
  optarglist_free( opts );

  Array *p = matrix_dtw_backtrack( m1 );
  
  opts =  optarglist( "weights=double*", weights);
  Array *avg = dtw_add_signals( a, b, C, opts);
  optarglist_free( opts );

  array_free( d );
  array_free( m1 );
  array_free( p );

  return avg;
}
Array* hierarchical_average_unequal_length( Array **data, const Array *distmat, 
                                                          SignalAverageFunctionUnequalLength avgfct, 
                                                          OptArgList *optargs ){
  Dendrogram *T, *Tsub; 
  ProgressBarFunction progress=NULL;  
  LinkageFunction linkage=dgram_dist_completelinkage;
  void *ptr;

  bool ismatrix;
  matrix_CHECKSQR( ismatrix, distmat );
  if( !ismatrix ) return NULL;


  optarg_PARSE_PTR( optargs, "linkage", linkage, LinkageFunction, ptr );
  optarg_PARSE_PTR( optargs, "progress", progress, ProgressBarFunction, ptr );

  int N,C;
  N = distmat->size[0];
  C = data[0]->size[1];

  dprintf("Data is (%i x %i)\n", N, C );

  /*------------------- computation --------------------------*/
  /* build hierarchical dendrogram */
  T = agglomerative_clustering( distmat, linkage );

  int i;
  double weights[2];
  int *indices = (int*)malloc( N*sizeof(int) );
  for( i=0; i<N; i++ ){
     indices[i] = 1;
  }  

  /* now walk the tree to find pairs of trials to match */
  Array *new, *avg;
  int trials_left = N;
  uint idx[2]={0,0};
  if( progress ) progress( PROGRESSBAR_INIT, N );


  Array **d = (Array*)malloc( N*sizeof(Array*) );
  for( i=0; i<N; i++ ){
     d[i] = data[i];
  }

  while( trials_left >= 2 ){
     if( progress ) progress( PROGRESSBAR_CONTINUE_LONG, N-trials_left );
     
     Tsub = dgram_get_deepest( T );
     idx[0] = Tsub->left->val;
     idx[1] = Tsub->right->val;
     
     weights[0] = indices[idx[0]]/(double)(indices[idx[0]]+indices[idx[1]]);
     weights[1] = indices[idx[1]]/(double)(indices[idx[0]]+indices[idx[1]]);

     new=avgfct( d, idx, weights, optargs );
     
     d[idx[0]]=new;
     d[idx[1]]=NULL;

     indices[idx[0]] += indices[idx[1]];
     indices[idx[1]]=-1;              /* never to be used again */
     
     /* replace node by leaf representing average(idx1, idx2) */
     Tsub->val = idx[0];
     Tsub->left = NULL;
     Tsub->right = NULL;    
     trials_left--;
  }
  /* idx[0] is the final average */
  avg = d[idx[0]];
  /*------------------- /computation --------------------------*/

  free( d );
  free( indices );

  return avg;
}

#endif