• API Main Page
  • Documentation
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

src/warping.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  *   Copyright (C) 2008 by Matthias Ihrke   *
00003  *   mihrke@uni-goettingen.de   *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019  ***************************************************************************/
00020 
00021 #include "warping.h"
00022 #include "optarg.h"
00023 #include "array.h"
00024 #include "linalg.h"
00025 #include "eeg.h"
00026 #include <math.h>
00027 
00034 double _slope_constraint( SlopeConstraint constraint, Array *d, Array *D, int i, int j ){
00035   double r;
00036 
00037   if( (i==0 || j==0) ){
00038      errprintf("i and j must be larger than zero '%i'\n", (int)constraint);
00039      return NAN;
00040   }
00041 
00042   r = mat_IDX( D, i-1, j-1 ) + 2*mat_IDX( d, i, j );
00043   switch( constraint ){
00044   case SLOPE_CONSTRAINT_NONE:
00045      r = MIN( r, mat_IDX( D, i, j-1) + mat_IDX( d, i, j ) );
00046      r = MIN( r, mat_IDX( D, i-1, j) + mat_IDX( d, i, j ) );
00047      break;
00048   case SLOPE_CONSTRAINT_LAX:
00049      if( j>=3 )
00050         r = MIN( r, mat_IDX( D, i-1, j-3 ) + 2*mat_IDX( d, i, j-2 ) + mat_IDX( d, i, j-1 ) + mat_IDX( d, i, j ) );
00051      if( j>=2 ){
00052         r = MIN( r, mat_IDX( D, i-1, j-2 ) + 2*mat_IDX( d, i, j-1 ) + mat_IDX( d, i, j ) );
00053      }
00054      if( i>=2 ){
00055         r = MIN( r, mat_IDX( D, i-2, j-1 ) + 2*mat_IDX( d, i-1, j ) + mat_IDX( d, i, j ) );
00056      }
00057      if( i>=3 ){
00058         r = MIN( r, mat_IDX( D, i-3, j-1 ) + 2*mat_IDX( d, i-2, j ) + mat_IDX( d, i-1, j ) + mat_IDX( d, i, j ) );
00059      }
00060      break;
00061   case SLOPE_CONSTRAINT_MEDIUM:
00062      if( j>=2 ){
00063         r = MIN( r, mat_IDX( D, i-1, j-2 ) + 2*mat_IDX( d, i, j-1 ) + mat_IDX( d, i, j ) );
00064      }
00065      if( i>=2 ){
00066         r = MIN( r, mat_IDX( D, i-2, j-1 ) + 2*mat_IDX( d, i-1, j ) + mat_IDX( d, i, j ) );
00067      }
00068      break;
00069   case SLOPE_CONSTRAINT_SEVERE:
00070      if( i>=2 && j>=3 ){
00071         r = MIN( r, mat_IDX( D, i-2, j-3 ) + 2*mat_IDX( d, i-1, j-2 ) + 2*mat_IDX( d, i, j-1 ) + mat_IDX( d, i, j ) );
00072      }
00073      if( j>=2 && i>=3 ){
00074         r = MIN( r, mat_IDX( D, i-3, j-2 ) + 2*mat_IDX( d, i-2, j-1 ) + 2*mat_IDX( d, i-1, j ) + mat_IDX( d, i, j ) );
00075      }
00076      break;
00077   default:
00078      errprintf("Slope constraint not supported '%i'\n", (int)constraint);
00079      r = NAN;
00080      break;
00081   }
00082   return r;
00083 }
00101 Array*      matrix_dtw_cumulate ( Array *mat, bool alloc, OptArgList *optargs ){
00102   /* j = 0,...,M-1
00103       k = 0,...,N-1
00104   */
00105   int j, k, M, N;
00106   double x;
00107   Array *D;
00108   SlopeConstraint constraint = SLOPE_CONSTRAINT_NONE;
00109 
00110   bool ismatrix;
00111   matrix_CHECK( ismatrix, mat );
00112   if( !ismatrix ) return NULL;
00113 
00114   if( optarglist_has_key( optargs, "slope_constraint" ) ){
00115      x = optarglist_scalar_by_key( optargs, "slope_constraint" );
00116      if( !isnan( x ) ) constraint=(SlopeConstraint)x;
00117   }
00118   dprintf("using slope constraint '%i'\n", constraint );
00119 
00120   D = array_copy( mat, TRUE );
00121 
00122   M = mat->size[0];
00123   N = mat->size[1];
00124    /* computing D_jk */
00125   for(j=0; j<M; j++){
00126     for(k=0; k<N; k++){
00127       if(k==0 && j==0) ;
00128       else if(k==0 && j>0){
00129           mat_IDX( D, j,k) = mat_IDX( D, j-1, k ) + mat_IDX( mat, j, k );
00130         } else if(k>0 && j==0){
00131           mat_IDX( D, j, k) = mat_IDX( D, j, k-1 ) + mat_IDX( mat, j, k );
00132         } else { /* j,k > 0 */
00133           mat_IDX( D, j, k) = _slope_constraint( constraint, mat, D, j, k ); 
00134         }
00135     }
00136   }
00137   if( !alloc ){
00138      memcpy( mat->data, D->data, D->nbytes );
00139      array_free( D );
00140      D = mat;
00141   }
00142      
00143   return D;
00144 }
00145 
00152 Array* matrix_dtw_backtrack ( const Array *d ){ 
00153   int i,j,k, M, N; /* j = 0,...,M-1
00154                           k = 0,...,N-1  */
00155   int idx;
00156   double left, down, downleft;
00157 
00158   bool isthing;
00159   matrix_CHECK( isthing, d );
00160   if( !isthing ) return NULL;
00161 
00162   M = d->size[0];
00163   N = d->size[1];
00164 
00165   Array *path,
00166      *path2 = array_new2( UINT, 2, 2, M+N ); /* dummy path */
00167 
00168   /*------------ computation -------------------------------*/
00169   /* Backtracking */
00170   j=M-1; k=N-1;
00171 
00172   idx = 1;
00173   array_INDEX2( path2, uint, 0, 0 ) = j;
00174   array_INDEX2( path2, uint, 1, 0 ) = k;
00175   while( j>0 || k>0 ){ 
00176      if( k==0 ){ /* base cases */
00177         j--;
00178      } else if( j==0 ){
00179         k--;
00180      } else { /* min( d[j-1][k], d[j-1][k-1], d[j][k-1] ) */
00181         left     = mat_IDX( d, j-1, k );
00182         down     = mat_IDX( d, j  , k-1 );
00183         downleft = mat_IDX( d, j-1, k-1 );
00184 
00185         (isnan( mat_IDX( d, j-1,k  ) ) )?(left=DBL_MAX):(left=mat_IDX( d, j-1, k ) );
00186         (isnan( mat_IDX( d, j  ,k-1) ) )?(down=DBL_MAX):(down=mat_IDX( d, j  ,k-1) );
00187         (isnan( mat_IDX( d, j-1,k-1) ) )?(downleft=DBL_MAX):(downleft=mat_IDX(d,j-1,k-1));
00188 
00189         /* this order of comparisons was chosen to ensure, that diagonal
00190             steps are preferred in case of equal distances
00191         */
00192         if( downleft<=left ){     /* downleft or down */
00193           if( downleft<=down ){
00194              k--; j--;                /* downleft */
00195           } else {
00196              k--;                         /* down */
00197           }
00198         } else {                          /* left or down */
00199           if( down<=left ){
00200              k--;                         /* down */
00201           } else {
00202              j--;                     /* left */
00203           }
00204         }
00205 
00206 
00207      }  /* if */
00208 
00209      array_INDEX2( path2, uint, 0, idx ) = j;
00210      array_INDEX2( path2, uint, 1, idx ) = k;
00211      idx++;
00212   }
00213   
00214   /* now the warppath has wrong order, reverse */ 
00215   path = array_new2( UINT, 2, 2, idx );
00216 
00217   for( i=0; i<idx; i++ ){
00218      array_INDEX2( path, uint, 0, idx-i-1 )=array_INDEX2( path2, uint, 0, i );
00219      array_INDEX2( path, uint, 1, idx-i-1 )=array_INDEX2( path2, uint, 1, i );
00220   }
00221   /*------------ /computation -------------------------------*/
00222 
00223   array_free( path2);
00224   return path;
00225 }
00226 
00242 Array*  dtw_add_signals( const Array *s1, const Array *s2, const Array *path, OptArgList *opts ){
00243   int ispath;
00244   warppath_CHECK( ispath, path );
00245   if( !ispath ) return NULL;
00246   
00247   if( s1->dtype!=DOUBLE || s2->dtype!=DOUBLE ){
00248      errprintf("Input signals must be double-arrays\n");
00249      return NULL;
00250   } 
00251   if( s1->ndim>2 || s2->ndim>2 || s1->ndim != s2->ndim ){
00252      warnprintf("Input signals should be 1D or 2D... continue at your own risk (%i,%i)\n",
00253                     s1->ndim, s2->ndim );
00254   }
00255   /* thin matrix-wrapper (in case of 1D data */
00256   Array *s1m, *s2m;
00257   s1m = array_fromptr2( DOUBLE, 2, s1->data, s1->size[0], (s1->ndim>1)?(s1->size[1]):1 );
00258   s2m = array_fromptr2( DOUBLE, 2, s2->data, s2->size[0], (s2->ndim>1)?(s2->size[1]):1 );
00259 
00260   if( s1m->size[1] != s2m->size[1] ){
00261      errprintf("Signals must have the same dimension 2\n");
00262      array_free( s1m );
00263      array_free( s2m );
00264      return NULL; 
00265   }
00266   int N1 = s1m->size[0];
00267   int N2 = s2m->size[0];
00268   int Nn = path->size[1];
00269   int p = s1m->size[1];
00270 
00271   double defaultw[2]={0.5,0.5};
00272   void *tmp;
00273   double *w=defaultw;
00274   optarg_PARSE_PTR( opts, "weights", w, double*, tmp );
00275   
00276 
00277   /* --------------------- computation ------------------------*/
00278   Array *wa = array_new2( DOUBLE, 2, Nn, p );
00279   int i,j;
00280   for( i=0; i<Nn; i++ ){
00281      for( j=0; j<p; j++ ){
00282         mat_IDX( wa, i, j )=w[0]*mat_IDX( s1m, array_INDEX2(path,uint,0,i),j )
00283           + w[1]*mat_IDX( s2m, array_INDEX2(path,uint,1,i),j );
00284      }
00285   }
00286   /* --------------------- /computation ------------------------*/
00287 
00288   array_free( s1m );
00289   array_free( s2m );
00290   return wa;
00291 }
00292 
00293 
00318 EEG* eeg_gibbons( EEG *eeg, int stimulus_marker, int response_marker, double k ){
00319 #ifdef FIXEEG
00320   double *times, *new_times;
00321   double meanRT=0, curRT;
00322   double step;
00323   int n, N, i, j, l;
00324   int numchan, c;
00325   int stim_onset=0, resp_onset=0;
00326   EEG *eeg_out;
00327 
00328   for( i=0; i<eeg->ntrials; i++ ){
00329      if( stimulus_marker>=eeg->nmarkers[i] ||
00330           response_marker>=eeg->nmarkers[i] ||
00331           stimulus_marker>=response_marker ){
00332         errprintf( "stimulus or response marker not correct (%i, %i) in trial %i, abort fct\n", 
00333                       stimulus_marker, response_marker, i );
00334         return NULL;
00335      }
00336   }
00337 
00338   eeg_out = eeg_init( eeg->nbchan, 1, eeg->n );
00339   eeg_out->sampling_rate=eeg->sampling_rate; 
00340   if( eeg->times ){
00341      eeg_out->times = (double*) malloc( eeg->n*sizeof(double) );
00342      memcpy( eeg_out->times, eeg->times, eeg->n*sizeof(double) );
00343   }
00344   if( eeg->chaninfo ){
00345      eeg_out->chaninfo = (ChannelInfo*) malloc( eeg->nbchan*sizeof(ChannelInfo) );
00346      memcpy( eeg_out->chaninfo, eeg->chaninfo,  eeg->nbchan*sizeof(ChannelInfo) );
00347   }
00348   eeg_append_comment( eeg_out, "output from eeg_gibbons()\n");
00349   eeg_init_markers( 2, eeg_out );
00350 
00351   /* for convenience */
00352   times = eeg->times;  
00353   N = eeg->ntrials;
00354   n = eeg->n;
00355   numchan = eeg->nbchan;
00356 
00357   /* gsl interpolation allocation */
00358   gsl_interp_accel *acc = gsl_interp_accel_alloc ();
00359   gsl_spline *spline = gsl_spline_alloc (gsl_interp_linear, n); /* gibbons uses linear */
00360 
00361   /* computing mean reaction time */
00362   for( i=0; i<N; i++ ){
00363      meanRT += times[eeg->markers[i][response_marker]] 
00364         - times[eeg->markers[i][stimulus_marker]];
00365   }
00366   meanRT /= (double)N; 
00367 
00368   /* compute warping for each trial */
00369   new_times = (double*) malloc( n*sizeof(double) );
00370   for( i=0; i<N; i++ ){
00371      memcpy( new_times, times, n*sizeof(double) ); /* copy times */
00372      stim_onset = eeg->markers[i][stimulus_marker];
00373      resp_onset = eeg->markers[i][response_marker];
00374      dprintf("markers = (%i,%i)\n", stim_onset, resp_onset);
00375      curRT = times[resp_onset]-times[stim_onset];
00376 
00377      for( j=stim_onset; j<resp_onset; j++ ){
00378         new_times[j] = times[j] + pow( times[j], k )
00379           /pow( curRT, k )*( meanRT-curRT ); /* gibbons eq. */
00380      }
00381      /* warp the rest of the segments linearly */
00382      step =(times[n-1]-curRT)/(double)(n-1-resp_onset);
00383      l = 0;
00384      for( j=resp_onset; j<n; j++ ){
00385         new_times[j]=meanRT + (l++)*step;
00386      }
00387 
00388      for( c=0; c<numchan; c++ ){ /* loop channels */
00389         gsl_spline_init (spline, new_times, eeg->data[c][i], n);
00390         for( j=0; j<n; j++ ){     /* and samples */
00391           eeg_out->data[c][0][j] += gsl_spline_eval ( spline, times[j], acc );
00392         }
00393      }
00394   }
00395   eeg_out->markers[0][0] = stim_onset;
00396   eeg_out->markers[0][1] = meanRT;
00397   eeg_out->marker_labels[0][0] = create_string( "stimulus" );
00398   eeg_out->marker_labels[0][1] = create_string( "response" );
00399 
00400   for( c=0; c<numchan; c++ ){
00401      for( j=0; j<n; j++ ){
00402         eeg_out->data[c][0][j] /= (double)N; /* average */
00403      }
00404   }
00405   
00406   /* free */
00407   gsl_spline_free (spline);
00408   gsl_interp_accel_free (acc);
00409   free( new_times );
00410 
00411   return eeg_out;
00412 #endif
00413 }
00414 
00415 /*****************************************************************************
00416 GOING TO BE OBSOLETE 
00417 ******************************************************************************/
00418 
00419 WarpPath* init_warppath(WarpPath *path, int n1, int n2){
00420   if(path==NULL){
00421      path = (WarpPath*)malloc(sizeof(WarpPath));
00422      path->t1=NULL;
00423      path->t2=NULL;
00424   }
00425   path->n = (n1+n2);
00426   path->n1= n1;
00427   path->n2= n2;
00428 
00429   if( path->t1 ){
00430      free( path->t1 );
00431   } 
00432   if( path->t2 ){
00433      free( path->t2 );
00434   }
00435   path->t1 = (int*)calloc(n1+n2, sizeof(int));
00436   path->t2 = (int*)calloc(n1+n2, sizeof(int));
00437 
00438   return path;
00439 }
00440 
00441 void      reset_warppath(WarpPath *P, int n1, int n2){
00442   memset( P->t1, 0, P->n*sizeof(int) );
00443   memset( P->t2, 0, P->n*sizeof(int) );
00444   P->n1 = n1; P->n2 = n2;
00445   P->n = 0;
00446   memset( P->t1, 0, n1*sizeof(int) );
00447   memset( P->t2, 0, n2*sizeof(int) );
00448 }
00449 
00450 void free_warppath(WarpPath *p){
00451     free(p->t1);
00452     free(p->t2);
00453     free(p);
00454 }
00455 
00456 void print_warppath( FILE *out, WarpPath *P ){
00457   int howmany, i;
00458   howmany=3;
00459 
00460   fprintf( out, "WarpPath '%p':\n"
00461               " n1 = %i\n"
00462               " n2 = %i\n"
00463               " n  = %i\n"
00464               " P  = ", P, P->n1, P->n2, P->n );
00465   if( P->n>=howmany ){
00466      for( i=0; i<howmany; i++ ){
00467         fprintf( out, " (%i -> %i),", P->t1[i], P->t2[i] );
00468      }
00469   } else if( P->n>0 ){
00470      for( i=0; i<P->n; i++ ){
00471         fprintf( out, " (%i -> %i),", P->t1[i], P->t2[i] );
00472      }
00473   } else {
00474      fprintf( out, "<empty>");
00475   }
00476   fprintf( out, "\n" );
00477 }
00478 
00488 void      dtw_cumulate_matrix  ( double **d, int M, int N, OptArgList *optargs ){
00489   /* j = 0,...,M-1
00490       k = 0,...,N-1
00491   */
00492   int j, k;
00493   double x;
00494   double **D;
00495   SlopeConstraint constraint = SLOPE_CONSTRAINT_NONE;
00496   
00497   if( optarglist_has_key( optargs, "slope_constraint" ) ){
00498      x = optarglist_scalar_by_key( optargs, "slope_constraint" );
00499      if( !isnan( x ) ) constraint=(SlopeConstraint)x;
00500   }
00501   dprintf("using slope constraint '%i'\n", constraint );
00502 
00503   D = dblpp_init( M, N );
00504   dblpp_copy( (const double**)d, D, M, N );
00505 
00506    /* computing D_jk */
00507   for(j=0; j<M; j++){
00508     for(k=0; k<N; k++){
00509       if(k==0 && j==0) ;
00510       else if(k==0 && j>0){
00511           D[j][k] = D[j-1][k] + d[j][k];
00512         } else if(k>0 && j==0){
00513           D[j][k] = D[j][k-1] + d[j][k];
00514         } else { /* j,k > 0 */
00515           D[j][k] = _slope_constraint( constraint, d, D, j, k ); // MIN(MIN(d[j][k-1], d[j-1][k]), d[j-1][k-1]);
00516         }
00517     }
00518   }
00519   
00520   dblpp_copy( (const double**)D, d, M, N );
00521   dblpp_free( D, M );
00522 }
00526 WarpPath* dtw_backtrack        ( const double **d, int M, int N, WarpPath *P ){ 
00527   int j,k; /* j = 0,...,M-1
00528                   k = 0,...,N-1  */
00529   int idx;
00530   double left, down, downleft;
00531 
00532   if( P==NULL ){
00533      P = init_warppath( NULL, M, N );
00534   }
00535 
00536   /* Backtracking */
00537   j=M-1; k=N-1;
00538 
00539   idx = 1;
00540   P->t1[0] = j;
00541   P->t2[0] = k;
00542   while( j>0 || k>0 ){
00543      if( k==0 ){ /* base cases */
00544         j--;
00545      } else if( j==0 ){
00546         k--;
00547      } else { /* min( d[j-1][k], d[j-1][k-1], d[j][k-1] ) */
00548       
00549         left     = d[j-1][k  ];
00550         down     = d[j  ][k-1];
00551         downleft = d[j-1][k-1];
00552 
00553         (isnan( d[j-1][k  ] ) )?(left     = DBL_MAX):(left     = d[j-1][k  ]);
00554         (isnan( d[j  ][k-1] ) )?(down     = DBL_MAX):(down     = d[j  ][k-1]);
00555         (isnan( d[j-1][k-1] ) )?(downleft = DBL_MAX):(downleft = d[j-1][k-1]);
00556 
00557 
00558         /* this order of comparisons was chosen to ensure, that diagonal
00559             steps are preferred in case of equal distances
00560         */
00561         if( downleft<=left ){     /* downleft or down */
00562           if( downleft<=down ){  
00563              k--; j--;                /* downleft */
00564           } else {
00565              k--;                         /* down */
00566           } 
00567         } else {                          /* left or down */
00568           if( down<=left ){
00569              k--;                         /* down */
00570           } else {
00571              j--;                     /* left */
00572           }
00573         }
00574 
00575 
00576      }  /* if */
00577 
00578      P->t1[idx] = j;
00579      P->t2[idx] = k;
00580      //  dprintf("(j, k), idx = (%i,%i), %i\n", j, k, idx);
00581      if(idx>0 && P->t1[idx]>=P->t1[idx-1] && P->t2[idx]>=P->t2[idx-1]  ){
00582         warnprintf("P=%p: idx=%i\n", P, idx);
00583      }
00584      idx++;
00585   }
00586   P->n=idx-1;
00587   //  print_warppath( stderr, P );
00588 
00589   /* now the warppath has wrong order, reverse */
00590   for( j=0; j<(P->n)/2; j++ ){
00591      swap2i( &(P->t1[j]), &(P->t1[P->n-1-j]) );
00592      swap2i( &(P->t2[j]), &(P->t2[P->n-1-j]) );
00593   }
00594   //  print_warppath( stderr, P );
00595 
00596   return P;
00597 }
00598 
00599 
00600 #ifdef EXPERIMENTAL
00601 
00606 double multidist( Array *dat, uint *idx ){  
00607   bool ismat;
00608   matrix_CHECK( ismat, dat );
00609   if( !ismat ) return NAN;
00610   
00611   double dist=0.0;
00612   int N=dat->size[0];
00613   int i,j;
00614 
00615   for( i=0; i<N-1; i++ ){
00616      for( j=i+1; j<N; j++ ){
00617         //      dprintf("(i,j)=%i,%i, idx=%i,%i\n",i,j, idx[i], idx[j]);
00618         dist += 0.9* ABS( mat_IDX( dat, i, idx[i] )-
00619                                 mat_IDX( dat, j, idx[j] ) ) +
00620           0.1* ABS( idx[i]-idx[j] );
00621      }
00622   }
00623   dist /= (double)((N*(N-1))/2);
00624 
00625   //  dprintf("d(i,j)=%f\n", dist );
00626   return dist;
00627 }
00628 
00631 Array* multiwarp( Array *in, uint window, OptArgList *optargs ){
00632   bool ismat;
00633   matrix_CHECK( ismat, in );
00634   if( !ismat ) return NULL;
00635   void *ptr;
00636   ProgressBarFunction progress=NULL;  
00637   optarg_PARSE_PTR( optargs, "progress", progress, ProgressBarFunction, ptr );
00638   if( progress )
00639      dprintf("found progressbar\n");
00640 
00641   int i, j;  
00642   int N=in->size[0]; /* num signals */
00643   int n=in->size[1]; /* num samples per signal */
00644 
00645   dprintf("N=%i,n=%i\n", N, n );
00646 
00647   /* allocate N-dimensional distance matrix */
00648   uint *dims;
00649   MALLOC( dims, N, uint );
00650   for( i=0; i<N; i++ )
00651      dims[i]=n;
00652   Array *distmat=array_new( DOUBLE, N, dims );
00653   free( dims );
00654 
00655   /* calculate distmat and cumulate at the same time */
00656   uint *idx, *cidx, *tidx;
00657   MALLOC( idx, N, uint );
00658   MALLOC( cidx, N, uint );
00659   MALLOC( tidx, N, uint );
00660 
00661   ulonglong lidx=0;
00662   Array *minvec=array_new2( DOUBLE, 1, (1<<N)-1 );
00663 
00664   if( progress )  progress( PROGRESSBAR_INIT, array_NUMEL( distmat ) );
00665   for( lidx=0; lidx<array_NUMEL( distmat ); lidx++ ){    
00666      /*--- distance calculation */
00667      array_calc_rowindex( lidx, distmat->size, distmat->ndim, idx );
00668      (*((double*)array_index( distmat, idx )))=multidist( in, idx );
00669      //  dprintf("lidx=%li, %i, %i, (1<<N)-1=%i\n", lidx, idx[0], idx[1], (1<<N)-1 );    
00670 
00671      /*--- cumulate the matrix */
00672      /* get all elements "surrounding" idx */
00673      for( i=0; i<( 1<<N )-1; i++ ){
00674         for( j=0; j<N; j++ ){ /* calculate index */
00675           cidx[j] = idx[j];
00676           if( CHECK_BIT( i+1, j ) ){
00677              cidx[j] = (cidx[j]==0)?(0):(cidx[j]-1);
00678           }
00679         }
00680         vec_IDX( minvec, i ) = (*(double*)array_index( distmat, cidx ));
00681      }
00682 
00683      /* cumulating */
00684      (*((double*)array_index( distmat, idx ))) += (*(double*)array_min( minvec ));
00685 
00686      if( progress && (lidx % n==0) & lidx/n % n==0 ) 
00687         progress( PROGRESSBAR_CONTINUE_LONG, lidx );
00688   }
00689   if( progress ) progress( PROGRESSBAR_FINISH, 0 );
00690 
00691   /*--- backtracking-phase */
00692   Array *p1=array_new2( UINT, 2, N, N*n ); /* dummy of maximum length */
00693 
00694   dprintf("Backtrack\n");
00695   lidx=0;
00696   for( i=0; i<N; i++ ){
00697      idx[i]=n-1; /* start in the upper right */
00698      (*(uint*)array_index2( p1, i, lidx ))=idx[i];
00699   }
00700   lidx=1;
00701   bool continue_flag, dontcheck_flag;
00702   double dmin;
00703   while(1){
00704      continue_flag=FALSE;
00705      dontcheck_flag=FALSE;
00706 
00707      /* get all elements "surrounding" idx and put minimal index in idx */
00708      dmin=DBL_MAX;
00709      for( i=0; i<( 1<<N )-1; i++ ){
00710         for( j=0; j<N; j++ ){ /* calculate index */
00711           cidx[j] = idx[j];
00712           if( CHECK_BIT( i+1, j ) ){
00713              if( cidx[j]==0 ) 
00714                 dontcheck_flag=TRUE;
00715              else 
00716                 cidx[j] = cidx[j]-1;
00717           }
00718         }
00719         if( !dontcheck_flag ){
00720           //          dprintf("Check field: (%i,%i,%i), dmin=%f\n", cidx[0], cidx[1], cidx[2], (*(double*)array_index( distmat, cidx )) );
00721 
00722           if(  (*(double*)array_index( distmat, cidx )) < dmin ){
00723              dmin = (*(double*)array_index( distmat, cidx ));
00724              for( j=0; j<N; j++ )
00725                 tidx[j]=cidx[j];
00726           }
00727         } else {
00728           dontcheck_flag=FALSE;
00729         }
00730      }   
00731      for( j=0; j<N; j++ )
00732         idx[j]=tidx[j];
00733 
00734      //  dprintf("Select field: (%i,%i,%i), dmin=%f\n", idx[0], idx[1], idx[2], dmin );
00735 
00736      /* advance path */
00737      for( j=0; j<N; j++ ){
00738         (*(uint*)array_index2( p1, j, lidx ))=idx[j];
00739      }
00740      lidx++;
00741 
00742      /* done? */
00743      for( i=0; i<N; i++ ){
00744         if( idx[i]>0 ){
00745           continue_flag=TRUE;
00746           break;
00747         }
00748      }
00749 
00750      if( !continue_flag )
00751         break;
00752   }
00753 
00754   /* reverse warppath and trim to correct size */
00755   Array *p2=array_new2( UINT, 2, N, lidx );
00756   for( i=0; i<lidx; i++ ){
00757      for( j=0; j<N; j++ ){
00758         array_INDEX2( p2, uint, j, lidx-i-1 ) = array_INDEX2( p1, uint, j, i );
00759      }
00760   }
00761 
00762   /* cleaning up */
00763   array_free( minvec );
00764   free( cidx );
00765   free( idx );
00766   free( tidx );
00767   array_free( distmat );
00768   array_free( p1 );
00769 
00770   return p2;
00771 }
00772 
00773 Array* multiwarp_add( Array *in, Array *times, Array *path ){
00774   bool ismat, isvec;
00775   matrix_CHECK( ismat, in );
00776   if( !ismat ) return NULL;
00777   vector_CHECK( isvec, times );
00778   if( !isvec ) return NULL;
00779 
00780   int i, j;  
00781   int N=in->size[0]; /* num signals */
00782   int n=in->size[1]; /* num samples per signal */
00783   int Nn=path->size[1];
00784 
00785   dprintf("N=%i,n=%i,Nn=%i\n", N, n, Nn );
00786 
00787   double *w;
00788   MALLOC( w, N, double );
00789   for( i=0; i<N; i++ ){
00790      w[i] = 1.0/(double)N;
00791   }
00792   dprintf("w=%f\n", w[0] );
00793 
00794   /* return times array in first dim */
00795   Array *wa = array_new2( DOUBLE, 2, 2, Nn );
00796   for( i=0; i<Nn; i++ ){
00797      for( j=0; j<N; j++ ){
00798         mat_IDX( wa, 0, i) += w[j]*vec_IDX( times, array_INDEX2(path,uint,j,i));
00799         mat_IDX( wa, 1, i) += w[j]*mat_IDX( in, j, array_INDEX2(path,uint,j,i));
00800      }
00801   }
00802   
00803   free( w );
00804   
00805   return wa;
00806 }
00807 
00808 
00809 #endif

Generated on Fri Jun 25 2010 14:10:21 for libeegtools by  doxygen 1.7.0