00001 /*************************************************************************** 00002 * Copyright (C) 2008/2009 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 "recurrence_plot.h" 00022 00051 Array* recplot( const Array *s1, const Array *s2, Array *out, double epsilon, 00052 OptArgList *optargs ){ 00053 int fan=-1; /* fixed amount of neighbours? */ 00054 double x; 00055 int i,j; 00056 if( optarglist_has_key( optargs, "fan" ) ){ 00057 x = optarglist_scalar_by_key( optargs, "fan" ); 00058 if( !isnan( x ) ) 00059 fan=(int)x; 00060 } 00061 if( !array_comparable( s1, s2 ) ){ 00062 return NULL; 00063 } 00064 if( s1->ndim>2 || s2->ndim>2 || s1->dtype!=DOUBLE || s2->dtype!=DOUBLE){ 00065 errprintf("Input arrays must be vectors or matrices\n"); 00066 return NULL; 00067 } 00068 00069 /* thin matrix-wrapper (in case of 1D data */ 00070 Array *s1m, *s2m; 00071 s1m = array_fromptr2( DOUBLE, 2, s1->data, s1->size[0], (s1->ndim>1)?(s1->size[1]):1 ); 00072 s2m = array_fromptr2( DOUBLE, 2, s2->data, s2->size[0], (s2->ndim>1)?(s2->size[1]):1 ); 00073 00074 int N = s1m->size[0]; 00075 int p = s1m->size[1]; 00076 if( !out ){ 00077 out = array_new2( DOUBLE, 2, N, N ); 00078 } else { 00079 if( out->ndim!=2 || out->dtype!=DOUBLE || out->size[0]<N || out->size[1]<N ){ 00080 errprintf("Output must be a N x N matrix (N=%i)\n", N); 00081 return out; 00082 } 00083 } 00084 00085 /* init FAN */ 00086 double *faneps=NULL; 00087 if( fan>0 ){ 00088 faneps=recplot_calculate_epsilons( s1m, s2m, NULL, fan ); 00089 } 00090 00091 /*---------- computation ----------- */ 00092 double eps; 00093 for( i=0; i<N; i++ ){ 00094 if( faneps ){ /* fixed amount of neighbours? */ 00095 eps = faneps[i]; 00096 } else { 00097 eps = epsilon; 00098 } 00099 00100 for( j=0; j<N; j++ ){ 00101 if( eps - vectordist_euclidean( (double*)array_INDEXMEM2( s1m, i, 0 ), 00102 (double*)array_INDEXMEM2( s2m, j, 0 ), 00103 p, optargs ) > 0 ){ 00104 mat_IDX(out, i, j) = 1.0; 00105 } else { 00106 mat_IDX(out, i, j) = 0.0; 00107 } 00108 } 00109 } 00110 /*---------- /computation ----------- */ 00111 00112 /* clean up */ 00113 array_free( s1m ); 00114 array_free( s2m ); 00115 if( faneps ) free( faneps ); 00116 return out; 00117 } 00118 00130 double* recplot_calculate_epsilons( Array *s1, Array *s2, double *eps, int fan ){ 00131 int i,j; 00132 bool ismatrix; 00133 if( !array_comparable( s1, s2 ) ){ 00134 return NULL; 00135 } 00136 matrix_CHECK( ismatrix, s1 ); 00137 if( !ismatrix ) return NULL; 00138 matrix_CHECK( ismatrix, s2 ); 00139 if( !ismatrix ) return NULL; 00140 00141 int N=s1->size[0]; 00142 int p=s1->size[1]; 00143 00144 if( !eps ){ 00145 eps = (double*)malloc( N*sizeof(double) ); 00146 } 00147 double *tmp = (double*)malloc( N*sizeof(double) ); 00148 double *epst= (double*)malloc( fan*sizeof(double) ); 00149 00150 /*---------- computation ----------- */ 00151 for( i=0; i<N; i++ ){ 00152 for( j=0; j<N; j++ ){ 00153 tmp[j] = vectordist_euclidean( (double*)array_INDEXMEM2( s1, i, j ), 00154 (double*)array_INDEXMEM2( s1, i, j ), 00155 p, NULL ); 00156 } 00157 /*int gsl_sort_smallest (double * dest, size_t k, const double * 00158 src, size_t stride, size_t n) 00159 This function copies the k smallest elements of the array src, 00160 of size n and stride stride, in ascending numerical order into 00161 the array dest. The size k of the subset must be less than or 00162 equal to n. The data src is not modified by this operation. */ 00163 gsl_sort_smallest( epst, fan, tmp, 1, N ); 00164 eps[i] = epst[fan-1]; 00165 } 00166 /*---------- /computation ----------- */ 00167 00168 free( tmp ); 00169 free( epst); 00170 return eps; 00171 } 00172 00173 /* /\** Calculate the Line-Of-Synchrony of a Cross-Recurrence-Plot, following */ 00174 /* the algorithm in Marwan et al. Cross recurrence plot based synchronization */ 00175 /* of time series. Nonlinear Processes in Geophysics (2002) vol. 9 (3-4) pp. 325-331. */ 00176 00177 /* \todo REPAIR THIS ALGORITHM! */ 00178 00179 /* Algorithm: */ 00180 /* \code */ 00181 /* Input: Parameters dx, dy */ 00182 /* 1. find recurrence point next to origin (i,j) */ 00183 /* 2. square window of size w=2 starting from (i,j); */ 00184 /* if another point in this window, goto 3. */ 00185 /* else w++, goto 2. */ 00186 /* 3. if point is found in x-direction, increase w in x-dir */ 00187 /* if point is found in y-direction, increase w in y-dir */ 00188 /* until: either no point is found or window is of */ 00189 /* size (w+dx, w+dy) */ 00190 /* after this step, window is of size (w+deltawx, w+deltawy) */ 00191 /* 4. set i^ = i + (w+deltawx)/2 */ 00192 /* set j^ = j + (w+deltawy)/2 */ 00193 /* 5. (i,j) = (i^,j^), goto 2. */ 00194 /* \endcode */ 00195 /* *\/ */ 00196 /* WarpPath* recplot_los_marwan( const RecurrencePlot *R, int dx, int dy ){ */ 00197 /* WarpPath *P,*P2; */ 00198 /* int ii, ij, /\* current LOS-point *\/ */ 00199 /* xcont,ycont,/\* continue in x/y direction flag *\/ */ 00200 /* deltawi, deltawj, /\* addition in x/y direction to (wi,wj) *\/ */ 00201 /* wi, wj; /\* dimensions of the window *\/ */ 00202 /* int i,j, pidx, breakflag; */ 00203 /* int no_recpoint_found, reached_border; */ 00204 00205 /* P = init_warppath( NULL, R->m, R->n ); */ 00206 00207 /* breakflag=0; */ 00208 /* ii=0; ij=0; */ 00209 /* for( i=0; MIN( R->m, R->n ); i++ ){ /\* find first point (1) *\/ */ 00210 /* for( j=0; j<=i; j++ ){ */ 00211 /* if( R->R[i][j] ){ */ 00212 /* ii = i; ij = j; */ 00213 /* breakflag=1; */ 00214 /* } else if( R->R[j][i] ){ */ 00215 /* ii = j; ij = i; */ 00216 /* breakflag=1; */ 00217 /* } */ 00218 /* if(breakflag) */ 00219 /* break; */ 00220 /* } */ 00221 /* if(breakflag) */ 00222 /* break; */ 00223 /* } */ 00224 /* /\* dprintf("(ii,ij) = (%i,%i)\n", ii,ij); *\/ */ 00225 /* pidx = 0; */ 00226 /* P->t1[pidx]=ii; */ 00227 /* P->t2[pidx]=ij; */ 00228 /* pidx++; */ 00229 00230 /* reached_border=0; */ 00231 /* while( !reached_border ){ */ 00232 /* wi = 2; */ 00233 /* wj = 2; */ 00234 /* no_recpoint_found=1; */ 00235 /* xcont = 0; */ 00236 /* ycont = 0; */ 00237 /* while( no_recpoint_found ){ /\* initial window sweep (2) *\/ */ 00238 /* if( ii+wi>=R->m || ij+wj>=R->n ){ */ 00239 /* reached_border=1; */ 00240 /* break; */ 00241 /* } */ 00242 00243 /* for( i=ii; i<ii+wi; i++ ){ */ 00244 /* if( R->R[i][ij+wj] ){ */ 00245 /* no_recpoint_found=0; */ 00246 /* xcont = 1; */ 00247 /* break; */ 00248 /* } */ 00249 /* } /\* for i *\/ */ 00250 00251 /* for( j=ij; j<ij+wj; j++ ){ */ 00252 /* if( R->R[ii+wi][j] ){ */ 00253 /* no_recpoint_found=0; */ 00254 /* ycont = 1; */ 00255 /* break; */ 00256 /* } */ 00257 /* } /\* for j *\/ */ 00258 00259 /* wi++; wj++; */ 00260 /* } /\* while(no_recpoint_found) *\/ */ 00261 /* if( reached_border ) */ 00262 /* break; */ 00263 00264 /* /\* dprintf("w=(%i,%i), cont=(%i,%i)\n", wi,wj, xcont, ycont ); *\/ */ 00265 /* deltawi=0; */ 00266 /* deltawj=0; */ 00267 /* while( xcont || ycont ){ /\* (3) *\/ */ 00268 /* if( ij+wj+deltawj>=(R->n-1) || ii+wi+deltawi>=(R->m-1) ){ */ 00269 /* reached_border=1; */ 00270 /* break; */ 00271 /* } */ 00272 /* if( xcont ){ /\* check border in x *\/ */ 00273 /* deltawi++; */ 00274 /* for( j=ij; j<ij+wj+deltawj; j++ ){ */ 00275 /* if( R->R[ii+wi+deltawi][j] ){ */ 00276 /* xcont = 1; */ 00277 /* } else { */ 00278 /* xcont = 0; */ 00279 /* } */ 00280 /* } /\* for j *\/ */ 00281 /* } */ 00282 00283 /* if( ycont ){ /\* check border in y *\/ */ 00284 /* deltawj++; */ 00285 /* for( i=ii; i<ii+wi+deltawi; i++ ){ */ 00286 /* if( R->R[i][ij+wj+deltawj] ){ */ 00287 /* ycont=1; */ 00288 /* } else { */ 00289 /* ycont=0; */ 00290 /* } */ 00291 /* } /\* for i *\/ */ 00292 /* } */ 00293 /* if( deltawi >= dx || deltawj >= dy ) */ 00294 /* break; */ 00295 /* } /\* while xcont *\/ */ 00296 /* /\* dprintf("deltaw=(%i,%i)\n", deltawi, deltawj ); *\/ */ 00297 /* if( reached_border ) */ 00298 /* break; */ 00299 00300 /* ii = ii + (wi+deltawi)/2; */ 00301 /* ij = ij + (wj+deltawj)/2; */ 00302 /* P->t1[pidx]=ii; */ 00303 /* P->t2[pidx]=ij; */ 00304 /* pidx++; */ 00305 /* } /\* while !reached_border *\/ */ 00306 /* P->t1[pidx]=R->m; */ 00307 /* P->t2[pidx]=R->n; */ 00308 /* P->n = ++pidx; */ 00309 00310 /* /\* fill gaps with bresenham *\/ */ 00311 /* P2 = init_warppath( NULL, R->m, R->n ); */ 00312 /* pidx = 0; */ 00313 /* int nump; */ 00314 /* int *line; */ 00315 /* line = (int*)malloc( R->m*R->n*2*sizeof(int)); /\* much too much *\/ */ 00316 00317 /* /\* for( i=0; i<P->n-1; i++ ){ *\/ */ 00318 /* /\* nump = bresenham_howmany_points( P->t1[i], P->t2[i], P->t1[i+1], P->t2[i+1] ); *\/ */ 00319 /* /\* line = bresenham( P->t1[i], P->t2[i], P->t1[i+1], P->t2[i+1], line ); *\/ */ 00320 /* /\* dprintf(" (%i,%i)->(%i,%i)\n", P->t1[i], P->t2[i], P->t1[i+1], P->t2[i+1] ); *\/ */ 00321 /* /\* for( j=0; j<nump; j++ ){ *\/ */ 00322 /* /\* dprintf(" line=(%i,%i)\n", line[(2*j)+0], line[(2*j)+1]); *\/ */ 00323 /* /\* P2->t1[pidx++]=line[(2*j)+0]; *\/ */ 00324 /* /\* P2->t2[pidx ]=line[(2*j)+1]; *\/ */ 00325 /* /\* } *\/ */ 00326 /* /\* } *\/ */ 00327 /* /\* P2->n = pidx; *\/ */ 00328 00329 /* free_warppath( P ); */ 00330 /* free( line ); */ 00331 00332 /* return P2; */ 00333 /* } */ 00334 00335 /* /\** Calculate the Line-Of-Synchrony of a Cross-Recurrence-Plot using */ 00336 /* a dynamic time-warping strategy. */ 00337 00338 /* Algorithm: */ 00339 /* \code */ 00340 /* Input: RecurrencePlot R (1 where recurrence, 0 otherwise) */ 00341 /* 1. calculate d = 1-R; */ 00342 /* 2. cumulate d, such that D[jk] = min{ D[j-1k], D[jk-1], D[j-1k-1] } */ 00343 /* 3. Backtrack */ 00344 /* \endcode */ 00345 /* (see \ref dtw for details) */ 00346 /* *\/ */ 00347 /* WarpPath* recplot_los_dtw( const RecurrencePlot *R ){ */ 00348 /* double **d; */ 00349 /* WarpPath *P; */ 00350 00351 00352 /* d = dblpp_init( R->m, R->n ); */ 00353 /* dblpp_copy( (const double**)R->R, d, R->m, R->n ); */ 00354 00355 /* /\* flip binary matrix *\/ */ 00356 /* scalar_minus_dblpp( 1.0, d, R->m, R->n ); */ 00357 /* /\* add a bias such that diagonal is preferred *\/ */ 00358 00359 00360 /* dtw_cumulate_matrix( d, R->m, R->n, NULL ); */ 00361 /* //dblpp_normalize_by_max( D, R->m, R->n ); */ 00362 /* P = dtw_backtrack( (const double**) d, R->m, R->n, NULL ); */ 00363 00364 /* dblpp_free( d, R->m ); */ 00365 00366 /* return P; */ 00367 /* } */ 00368 00369 /* /\** Calculate the Line-Of-Synchrony of a Cross-Recurrence-Plot using */ 00370 /* a regularized dynamic time-warping strategy. */ 00371 00372 /* Algorithm: */ 00373 /* \code */ 00374 /* Input: RecurrencePlot R (1 where recurrence, 0 otherwise) */ 00375 /* 1. calculate d = 1-R; */ 00376 /* 2. get the regularization function G that is the distance transform of the */ 00377 /* linear interpolation between corresponding points in time */ 00378 /* 3. d = d*G; */ 00379 /* 4. cumulate d, such that D[jk] = min{ D[j-1k], D[jk-1], D[j-1k-1] } */ 00380 /* 5. Backtrack */ 00381 /* \endcode */ 00382 /* (see \ref dtw for details) */ 00383 00384 /* \param R the recurrence plot */ 00385 /* \param markers 2xnmarkers array; markers[0] is markers of 1st signal, */ 00386 /* markers[1] for second; if NULL, { {0,n-1}, {0,n-1} } is used */ 00387 /* \param nmarkers number of markers */ 00388 /* \return the warping function */ 00389 /* *\/ */ 00390 /* WarpPath* recplot_los_dtw_markers( const RecurrencePlot *R, int **markers, int nmarkers ){ */ 00391 /* double **d, **G; */ 00392 /* WarpPath *P; */ 00393 /* int mflag=0; */ 00394 /* int n; */ 00395 00396 /* n = R->m; */ 00397 /* d = dblpp_init( R->m, R->n ); */ 00398 /* dblpp_copy( (const double**)R->R, d, R->m, R->n ); */ 00399 00400 /* /\* flip binary matrix *\/ */ 00401 /* scalar_minus_dblpp( 1.0, d, R->m, R->n ); */ 00402 00403 /* /\* add a bias such that f is preferred *\/ */ 00404 /* if( !markers ){ */ 00405 /* mflag = 1; */ 00406 /* markers = (int**)malloc( 2*sizeof( int* ) ); */ 00407 /* markers[0] = (int*)malloc( 2*sizeof( int ) ); */ 00408 /* markers[1] = (int*)malloc( 2*sizeof( int ) ); */ 00409 /* markers[0][0] = 0; markers[1][0] = 0; */ 00410 /* markers[0][1] = n-1; markers[1][1] = n-1; */ 00411 /* } */ 00412 /* OptArgList *opts = optarglist( "markers=int**,nmarkers=int", */ 00413 /* markers, nmarkers ); */ 00414 /* G = regularization_linear_points( ALLOC_IN_FCT, R->m, opts ); */ 00415 /* optarglist_free( opts ); */ 00416 /* dblpp_normalize_by_max( G, R->m, R->n ); */ 00417 00418 /* /\* apply *\/ */ 00419 /* dblpp_add_dblpp( d, (const double**)G, R->m, R->n ); */ 00420 00421 /* /\* dtw *\/ */ 00422 /* dtw_cumulate_matrix( d, R->m, R->n, NULL ); */ 00423 /* P = dtw_backtrack( (const double**) d, R->m, R->n, NULL ); */ 00424 00425 /* /\* clean *\/ */ 00426 /* dblpp_free( d, R->m ); */ 00427 /* dblpp_free( G, R->m ); */ 00428 /* if( mflag ){ */ 00429 /* free( markers[0] ); */ 00430 /* free( markers[1] ); */ 00431 /* free( markers ); */ 00432 /* } */ 00433 00434 /* return P; */ 00435 /* } */ 00436 00437 /* /\** Calculate the Line-Of-Synchrony of a Cross-Recurrence-Plot using */ 00438 /* an algorithm based on the distance transform of the plot. */ 00439 00440 /* Algorithm: */ 00441 /* \code */ 00442 /* Input: RecurrencePlot R (1 where recurrence, 0 otherwise) */ 00443 /* 1. calculate d = DT(R) \ref disttransform */ 00444 /* 2. cumulate d, such that D[jk] = min{ D[j-1k], D[jk-1], D[j-1k-1] } */ 00445 /* 3. Backtrack */ 00446 /* \endcode */ 00447 /* (see \ref dtw for details) */ 00448 00449 /* \param R the recurrence plot */ 00450 /* \return the warping function */ 00451 /* *\/ */ 00452 /* WarpPath* recplot_los_disttransform( const RecurrencePlot *R ){ */ 00453 /* WarpPath *P; */ 00454 /* int **I; */ 00455 /* double **d; */ 00456 /* int i,j; */ 00457 00458 /* /\* dummy, needed for distance transform *\/ */ 00459 /* I = (int**) malloc( R->m*sizeof( int* ) ); */ 00460 /* for( i=0; i<R->m; i++ ){ */ 00461 /* I[i] = (int*) malloc( R->n*sizeof( int ) ); */ 00462 /* for( j=0; j<R->n; j++ ){ */ 00463 /* I[i][j] = 0; */ 00464 /* if( R->R[i][j]>0 ){ */ 00465 /* I[i][j] = 1; */ 00466 /* } */ 00467 /* } */ 00468 /* } */ 00469 00470 /* d = disttransform_deadreckoning( I, R->m, R->n, ALLOC_IN_FCT ); */ 00471 /* dtw_cumulate_matrix( d, R->m, R->n, NULL ); */ 00472 /* P = dtw_backtrack( (const double**) d, R->m, R->n, NULL ); */ 00473 00474 /* /\* free *\/ */ 00475 /* dblpp_free( d, R->m ); */ 00476 /* for( i=0; i<R->m; i++ ) */ 00477 /* free( I[i] ); */ 00478 /* free( I ); */ 00479 00480 /* return P; */ 00481 /* } */ 00482 00483 /* /\** Calculate the Line-Of-Synchrony of a Cross-Recurrence-Plot using */ 00484 /* a dynamic time-warping strategy. */ 00485 00486 /* Algorithm: */ 00487 /* \code */ 00488 /* Input: RecurrencePlot R (1 where recurrence, 0 otherwise) */ 00489 /* 1. calculate d = 2-R; */ 00490 /* 2. add a small amount of noise: */ 00491 /* if d_jk = 1, d_jk = d_jk + epsilon */ 00492 /* 2. cumulate d, such that D[jk] = min{ D[j-1k], D[jk-1], D[j-1k-1] } */ 00493 /* 3. Backtrack */ 00494 /* \endcode */ 00495 /* (see \ref dtw for details) */ 00496 /* *\/ */ 00497 /* WarpPath* recplot_los_dtw_noise( const RecurrencePlot *R ){ */ 00498 /* double **d; */ 00499 /* WarpPath *P; */ 00500 /* double noise; */ 00501 /* double noiseamp = 0.01; */ 00502 /* int i,j; */ 00503 00504 /* srand((long)time(NULL)); */ 00505 /* d = dblpp_init( R->m, R->n ); */ 00506 /* dblpp_copy( (const double**)R->R, d, R->m, R->n ); */ 00507 00508 /* /\* flip binary matrix *\/ */ 00509 /* scalar_minus_dblpp( 1.0, d, R->m, R->n ); */ 00510 00511 /* /\* add small amount of noise *\/ */ 00512 /* for( i=0; i<R->m; i++ ){ */ 00513 /* for( j=0; j<R->n; j++ ){ */ 00514 /* noise = noiseamp*(((double)rand()) / RAND_MAX); */ 00515 /* // if( d[i][j]<2 ) */ 00516 /* d[i][j] += noise; */ 00517 /* } */ 00518 /* } */ 00519 00520 /* OptArgList *opt = optarglist( "slope_constraint=int", SLOPE_CONSTRAINT_LAX ); */ 00521 /* dtw_cumulate_matrix( d, R->m, R->n, opt ); */ 00522 /* optarglist_free( opt ); */ 00523 00524 /* P = dtw_backtrack( (const double**) d, R->m, R->n, NULL ); */ 00525 00526 /* dblpp_free( d, R->m ); */ 00527 00528 /* return P; */ 00529 /* } */
1.7.0