1 #include "mrilib.h"
2
3 typedef struct {
4 int bot , top ;
5 float pval ;
6 } spanfit ;
7
8 float evaluate_span( int , int , int , int , float * , float ** ) ;
9
10 /*-------------------------------------------------------------------*/
11
find_best_span(int ndim,int nvec,int minspan,float * cvec,float ** bvec)12 spanfit find_best_span( int ndim , int nvec , int minspan ,
13 float * cvec , float ** bvec )
14 {
15 spanfit result = {0,0,0.0} ;
16 int ii,kk , bot,top , bot_best,top_best ;
17 float val , val_best ;
18
19 if( minspan < 3 || ndim < minspan || nvec < 100 ) return result ;
20 if( cvec == NULL || bvec == NULL ) return result ;
21
22 val_best = -1.0 ;
23 for( bot=0 ; bot < ndim+1-minspan ; bot++ ){
24 for( top=bot+minspan-1 ; top < ndim ; top++ ){
25 val = evaluate_span( ndim,nvec , bot,top , cvec,bvec ) ;
26 if( val > val_best ){
27 val_best = val ; bot_best = bot ; top_best = top ;
28 }
29 }
30 }
31
32 evaluate_span( 0,0,0,0,NULL,NULL ) ;
33
34 result.bot = bot_best ; result.top = top_best ; result.pval = val_best ;
35 return result ;
36 }
37
38 /*-------------------------------------------------------------------*/
39
evaluate_span(int ndim,int nvec,int bot,int top,float * cvec,float ** bvec)40 float evaluate_span( int ndim, int nvec,
41 int bot , int top , float * cvec , float ** bvec )
42 {
43 int kk , ibot=bot,itop=top , nneg ;
44 float cbar, *qvec ;
45 register int ii ;
46 register float sum ;
47
48 static float * bsum=NULL , * cnorm=NULL ;
49 static int nbsum=0 , ncnorm=0 ;
50
51 if( nvec > nbsum ){
52 if( bsum != NULL ) free(bsum) ;
53 bsum = (float *) malloc(sizeof(float)*nvec) ;
54 nbsum = nvec ;
55 } else if( nvec <= 0 ){
56 if( bsum != NULL ){ free(bsum) ; bsum = NULL; nbsum = 0; }
57 if( cnorm != NULL ){ free(cnorm); cnorm = NULL; ncnorm = 0; }
58 return 0.0 ;
59 }
60
61 if( ndim > ncnorm ){
62 if( cnorm != NULL ) free(cnorm) ;
63 cnorm = (float *) malloc(sizeof(float)*ndim) ;
64 ncnorm = ndim ;
65 }
66
67 /* compute cnorm = cvec-cbar */
68
69 sum = 0.0 ;
70 for( ii=ibot ; ii <= itop ; ii++ ) sum += cvec[ii] ;
71 cbar = sum/(itop-ibot+1) ; sum = 0.0 ;
72 for( ii=ibot ; ii <= itop ; ii++ ){
73 cnorm[ii] = cvec[ii] - cbar ; sum += cnorm[ii]*cnorm[ii] ;
74 }
75 if( sum <= 0.0 ) return 1.0 ; /* [cvec-cbar]=0 is perfect */
76
77 /* project each bvec onto cnorm */
78
79 for( kk=0 ; kk < nvec ; kk++ ){
80 qvec = bvec[kk] ; sum = 0.0 ;
81 for( ii=ibot ; ii <= itop ; ii++ ) sum += qvec[ii] * cnorm[ii] ;
82 bsum[kk] = sum ;
83 }
84
85 /* find number of bsums less than 0 */
86
87 for( nneg=ii=0 ; ii < nvec ; ii++ ) if( bsum[ii] <= 0.0 ) nneg++ ;
88
89 /* return value is fraction of negative bsum values */
90
91 sum = (float)nneg / (float)nvec ;
92 return sum ;
93 }
94