00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
00103
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
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 {
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;
00154
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 );
00167
00168
00169
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 ){
00177 j--;
00178 } else if( j==0 ){
00179 k--;
00180 } else {
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
00190
00191
00192 if( downleft<=left ){
00193 if( downleft<=down ){
00194 k--; j--;
00195 } else {
00196 k--;
00197 }
00198 } else {
00199 if( down<=left ){
00200 k--;
00201 } else {
00202 j--;
00203 }
00204 }
00205
00206
00207 }
00208
00209 array_INDEX2( path2, uint, 0, idx ) = j;
00210 array_INDEX2( path2, uint, 1, idx ) = k;
00211 idx++;
00212 }
00213
00214
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
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
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
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
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
00352 times = eeg->times;
00353 N = eeg->ntrials;
00354 n = eeg->n;
00355 numchan = eeg->nbchan;
00356
00357
00358 gsl_interp_accel *acc = gsl_interp_accel_alloc ();
00359 gsl_spline *spline = gsl_spline_alloc (gsl_interp_linear, n);
00360
00361
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
00369 new_times = (double*) malloc( n*sizeof(double) );
00370 for( i=0; i<N; i++ ){
00371 memcpy( new_times, times, n*sizeof(double) );
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 );
00380 }
00381
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++ ){
00389 gsl_spline_init (spline, new_times, eeg->data[c][i], n);
00390 for( j=0; j<n; j++ ){
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;
00403 }
00404 }
00405
00406
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
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
00490
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
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 {
00515 D[j][k] = _slope_constraint( constraint, d, D, j, k );
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;
00528
00529 int idx;
00530 double left, down, downleft;
00531
00532 if( P==NULL ){
00533 P = init_warppath( NULL, M, N );
00534 }
00535
00536
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 ){
00544 j--;
00545 } else if( j==0 ){
00546 k--;
00547 } else {
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
00559
00560
00561 if( downleft<=left ){
00562 if( downleft<=down ){
00563 k--; j--;
00564 } else {
00565 k--;
00566 }
00567 } else {
00568 if( down<=left ){
00569 k--;
00570 } else {
00571 j--;
00572 }
00573 }
00574
00575
00576 }
00577
00578 P->t1[idx] = j;
00579 P->t2[idx] = k;
00580
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
00588
00589
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
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
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
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];
00643 int n=in->size[1];
00644
00645 dprintf("N=%i,n=%i\n", N, n );
00646
00647
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
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
00667 array_calc_rowindex( lidx, distmat->size, distmat->ndim, idx );
00668 (*((double*)array_index( distmat, idx )))=multidist( in, idx );
00669
00670
00671
00672
00673 for( i=0; i<( 1<<N )-1; i++ ){
00674 for( j=0; j<N; j++ ){
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
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
00692 Array *p1=array_new2( UINT, 2, N, N*n );
00693
00694 dprintf("Backtrack\n");
00695 lidx=0;
00696 for( i=0; i<N; i++ ){
00697 idx[i]=n-1;
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
00708 dmin=DBL_MAX;
00709 for( i=0; i<( 1<<N )-1; i++ ){
00710 for( j=0; j<N; j++ ){
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
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
00735
00736
00737 for( j=0; j<N; j++ ){
00738 (*(uint*)array_index2( p1, j, lidx ))=idx[j];
00739 }
00740 lidx++;
00741
00742
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
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
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];
00782 int n=in->size[1];
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
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