| 1 | #include "mltaln.h" |
|---|
| 2 | #include "dp.h" |
|---|
| 3 | |
|---|
| 4 | #define DEBUG 0 |
|---|
| 5 | |
|---|
| 6 | static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize ) |
|---|
| 7 | { |
|---|
| 8 | int j, k, l; |
|---|
| 9 | float scarr[26]; |
|---|
| 10 | float **cpmxpd = floatwork; |
|---|
| 11 | int **cpmxpdn = intwork; |
|---|
| 12 | int count = 0; |
|---|
| 13 | |
|---|
| 14 | if( initialize ) |
|---|
| 15 | { |
|---|
| 16 | for( j=0; j<lgth2; j++ ) |
|---|
| 17 | { |
|---|
| 18 | count = 0; |
|---|
| 19 | for( l=0; l<26; l++ ) |
|---|
| 20 | { |
|---|
| 21 | if( cpmx2[l][j] ) |
|---|
| 22 | { |
|---|
| 23 | cpmxpd[count][j] = cpmx2[l][j]; |
|---|
| 24 | cpmxpdn[count][j] = l; |
|---|
| 25 | count++; |
|---|
| 26 | } |
|---|
| 27 | } |
|---|
| 28 | cpmxpdn[count][j] = -1; |
|---|
| 29 | } |
|---|
| 30 | } |
|---|
| 31 | |
|---|
| 32 | for( l=0; l<26; l++ ) |
|---|
| 33 | { |
|---|
| 34 | scarr[l] = 0.0; |
|---|
| 35 | for( k=0; k<26; k++ ) |
|---|
| 36 | scarr[l] += n_dis[k][l] * cpmx1[k][i1]; |
|---|
| 37 | } |
|---|
| 38 | for( j=0; j<lgth2; j++ ) |
|---|
| 39 | { |
|---|
| 40 | match[j] = 0; |
|---|
| 41 | for( k=0; cpmxpdn[k][j] > -1; k++ ) |
|---|
| 42 | match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j]; |
|---|
| 43 | } |
|---|
| 44 | } |
|---|
| 45 | |
|---|
| 46 | static float Atracking( float *lasthorizontalw, float *lastverticalw, |
|---|
| 47 | char **seq1, char **seq2, |
|---|
| 48 | char **mseq1, char **mseq2, |
|---|
| 49 | float **cpmx1, float **cpmx2, |
|---|
| 50 | int **ijp, int icyc, int jcyc ) |
|---|
| 51 | { |
|---|
| 52 | int i, j, k, l, iin, jin, ifi, jfi, lgth1, lgth2; |
|---|
| 53 | // char gap[] = "-"; |
|---|
| 54 | char *gap; |
|---|
| 55 | float wm; |
|---|
| 56 | gap = newgapstr; |
|---|
| 57 | lgth1 = strlen( seq1[0] ); |
|---|
| 58 | lgth2 = strlen( seq2[0] ); |
|---|
| 59 | |
|---|
| 60 | #if DEBUG |
|---|
| 61 | for( i=0; i<lgth1; i++ ) |
|---|
| 62 | { |
|---|
| 63 | fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] ); |
|---|
| 64 | } |
|---|
| 65 | #endif |
|---|
| 66 | |
|---|
| 67 | if( outgap == 1 ) |
|---|
| 68 | ; |
|---|
| 69 | else |
|---|
| 70 | { |
|---|
| 71 | wm = lastverticalw[0]; |
|---|
| 72 | for( i=0; i<lgth1; i++ ) |
|---|
| 73 | { |
|---|
| 74 | if( lastverticalw[i] >= wm ) |
|---|
| 75 | { |
|---|
| 76 | wm = lastverticalw[i]; |
|---|
| 77 | iin = i; jin = lgth2-1; |
|---|
| 78 | ijp[lgth1][lgth2] = +( lgth1 - i ); |
|---|
| 79 | } |
|---|
| 80 | } |
|---|
| 81 | for( j=0; j<lgth2; j++ ) |
|---|
| 82 | { |
|---|
| 83 | if( lasthorizontalw[j] >= wm ) |
|---|
| 84 | { |
|---|
| 85 | wm = lasthorizontalw[j]; |
|---|
| 86 | iin = lgth1-1; jin = j; |
|---|
| 87 | ijp[lgth1][lgth2] = -( lgth2 - j ); |
|---|
| 88 | } |
|---|
| 89 | } |
|---|
| 90 | } |
|---|
| 91 | |
|---|
| 92 | for( i=0; i<lgth1+1; i++ ) |
|---|
| 93 | { |
|---|
| 94 | ijp[i][0] = i + 1; |
|---|
| 95 | } |
|---|
| 96 | for( j=0; j<lgth2+1; j++ ) |
|---|
| 97 | { |
|---|
| 98 | ijp[0][j] = -( j + 1 ); |
|---|
| 99 | } |
|---|
| 100 | |
|---|
| 101 | for( i=0; i<icyc; i++ ) |
|---|
| 102 | { |
|---|
| 103 | mseq1[i] += lgth1+lgth2; |
|---|
| 104 | *mseq1[i] = 0; |
|---|
| 105 | } |
|---|
| 106 | for( j=0; j<jcyc; j++ ) |
|---|
| 107 | { |
|---|
| 108 | mseq2[j] += lgth1+lgth2; |
|---|
| 109 | *mseq2[j] = 0; |
|---|
| 110 | } |
|---|
| 111 | iin = lgth1; jin = lgth2; |
|---|
| 112 | for( k=0; k<=lgth1+lgth2; k++ ) |
|---|
| 113 | { |
|---|
| 114 | if( ijp[iin][jin] < 0 ) |
|---|
| 115 | { |
|---|
| 116 | ifi = iin-1; jfi = jin+ijp[iin][jin]; |
|---|
| 117 | } |
|---|
| 118 | else if( ijp[iin][jin] > 0 ) |
|---|
| 119 | { |
|---|
| 120 | ifi = iin-ijp[iin][jin]; jfi = jin-1; |
|---|
| 121 | } |
|---|
| 122 | else |
|---|
| 123 | { |
|---|
| 124 | ifi = iin-1; jfi = jin-1; |
|---|
| 125 | } |
|---|
| 126 | l = iin - ifi; |
|---|
| 127 | while( --l ) |
|---|
| 128 | { |
|---|
| 129 | for( i=0; i<icyc; i++ ) |
|---|
| 130 | *--mseq1[i] = seq1[i][ifi+l]; |
|---|
| 131 | for( j=0; j<jcyc; j++ ) |
|---|
| 132 | *--mseq2[j] = *gap; |
|---|
| 133 | k++; |
|---|
| 134 | } |
|---|
| 135 | l= jin - jfi; |
|---|
| 136 | while( --l ) |
|---|
| 137 | { |
|---|
| 138 | for( i=0; i<icyc; i++ ) |
|---|
| 139 | *--mseq1[i] = *gap; |
|---|
| 140 | for( j=0; j<jcyc; j++ ) |
|---|
| 141 | *--mseq2[j] = seq2[j][jfi+l]; |
|---|
| 142 | k++; |
|---|
| 143 | } |
|---|
| 144 | if( iin <= 0 || jin <= 0 ) break; |
|---|
| 145 | for( i=0; i<icyc; i++ ) |
|---|
| 146 | *--mseq1[i] = seq1[i][ifi]; |
|---|
| 147 | for( j=0; j<jcyc; j++ ) |
|---|
| 148 | *--mseq2[j] = seq2[j][jfi]; |
|---|
| 149 | k++; |
|---|
| 150 | iin = ifi; jin = jfi; |
|---|
| 151 | } |
|---|
| 152 | return( 0.0 ); |
|---|
| 153 | } |
|---|
| 154 | |
|---|
| 155 | |
|---|
| 156 | float Aalign( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen ) |
|---|
| 157 | /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */ |
|---|
| 158 | { |
|---|
| 159 | register int i, j; |
|---|
| 160 | int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */ |
|---|
| 161 | int lgth1, lgth2; |
|---|
| 162 | int resultlen; |
|---|
| 163 | float wm = 0.0; /* int ?????? */ |
|---|
| 164 | float g; |
|---|
| 165 | float x; |
|---|
| 166 | static TLS float mi, *m; |
|---|
| 167 | static TLS int **ijp; |
|---|
| 168 | static TLS int mpi, *mp; |
|---|
| 169 | static TLS float *currentw; |
|---|
| 170 | static TLS float *previousw; |
|---|
| 171 | static TLS float *match; |
|---|
| 172 | static TLS float *initverticalw; /* kufuu sureba iranai */ |
|---|
| 173 | static TLS float *lastverticalw; /* kufuu sureba iranai */ |
|---|
| 174 | static TLS char **mseq1; |
|---|
| 175 | static TLS char **mseq2; |
|---|
| 176 | static TLS char **mseq; |
|---|
| 177 | static TLS float **cpmx1; |
|---|
| 178 | static TLS float **cpmx2; |
|---|
| 179 | static TLS int **intwork; |
|---|
| 180 | static TLS float **floatwork; |
|---|
| 181 | static TLS int orlgth1 = 0, orlgth2 = 0; |
|---|
| 182 | |
|---|
| 183 | #if DEBUG |
|---|
| 184 | fprintf( stderr, "eff in SA+++align\n" ); |
|---|
| 185 | for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] ); |
|---|
| 186 | #endif |
|---|
| 187 | if( orlgth1 == 0 ) |
|---|
| 188 | { |
|---|
| 189 | mseq1 = AllocateCharMtx( njob, 1 ); |
|---|
| 190 | mseq2 = AllocateCharMtx( njob, 1 ); /* by J. Thompson */ |
|---|
| 191 | } |
|---|
| 192 | |
|---|
| 193 | lgth1 = strlen( seq1[0] ); |
|---|
| 194 | lgth2 = strlen( seq2[0] ); |
|---|
| 195 | |
|---|
| 196 | if( lgth1 > orlgth1 || lgth2 > orlgth2 ) |
|---|
| 197 | { |
|---|
| 198 | int ll1, ll2; |
|---|
| 199 | |
|---|
| 200 | if( orlgth1 > 0 && orlgth2 > 0 ) |
|---|
| 201 | { |
|---|
| 202 | FreeFloatVec( currentw ); |
|---|
| 203 | FreeFloatVec( previousw ); |
|---|
| 204 | FreeFloatVec( match ); |
|---|
| 205 | FreeFloatVec( initverticalw ); |
|---|
| 206 | FreeFloatVec( lastverticalw ); |
|---|
| 207 | |
|---|
| 208 | FreeFloatVec( m ); |
|---|
| 209 | FreeIntVec( mp ); |
|---|
| 210 | |
|---|
| 211 | FreeCharMtx( mseq ); |
|---|
| 212 | |
|---|
| 213 | FreeFloatMtx( cpmx1 ); |
|---|
| 214 | FreeFloatMtx( cpmx2 ); |
|---|
| 215 | |
|---|
| 216 | FreeFloatMtx( floatwork ); |
|---|
| 217 | FreeIntMtx( intwork ); |
|---|
| 218 | } |
|---|
| 219 | |
|---|
| 220 | ll1 = MAX( (int)(1.1*lgth1), orlgth1 ) + 100; |
|---|
| 221 | ll2 = MAX( (int)(1.1*lgth2), orlgth2 ) + 100; |
|---|
| 222 | |
|---|
| 223 | fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 ); |
|---|
| 224 | |
|---|
| 225 | currentw = AllocateFloatVec( ll2+2 ); |
|---|
| 226 | previousw = AllocateFloatVec( ll2+2 ); |
|---|
| 227 | match = AllocateFloatVec( ll2+2 ); |
|---|
| 228 | |
|---|
| 229 | initverticalw = AllocateFloatVec( ll1+2 ); |
|---|
| 230 | lastverticalw = AllocateFloatVec( ll1+2 ); |
|---|
| 231 | |
|---|
| 232 | m = AllocateFloatVec( ll2+2 ); |
|---|
| 233 | mp = AllocateIntVec( ll2+2 ); |
|---|
| 234 | |
|---|
| 235 | mseq = AllocateCharMtx( njob, ll1+ll2 ); |
|---|
| 236 | |
|---|
| 237 | cpmx1 = AllocateFloatMtx( 26, ll1+2 ); |
|---|
| 238 | cpmx2 = AllocateFloatMtx( 26, ll2+2 ); |
|---|
| 239 | |
|---|
| 240 | floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); |
|---|
| 241 | intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); |
|---|
| 242 | |
|---|
| 243 | fprintf( stderr, "succeeded\n" ); |
|---|
| 244 | |
|---|
| 245 | orlgth1 = ll1; |
|---|
| 246 | orlgth2 = ll2; |
|---|
| 247 | } |
|---|
| 248 | |
|---|
| 249 | for( i=0; i<icyc; i++ ) mseq1[i] = mseq[i]; |
|---|
| 250 | for( j=0; j<jcyc; j++ ) mseq2[j] = mseq[icyc+j]; |
|---|
| 251 | |
|---|
| 252 | |
|---|
| 253 | if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 ) |
|---|
| 254 | { |
|---|
| 255 | int ll1, ll2; |
|---|
| 256 | |
|---|
| 257 | if( commonAlloc1 && commonAlloc2 ) |
|---|
| 258 | { |
|---|
| 259 | FreeIntMtx( commonIP ); |
|---|
| 260 | } |
|---|
| 261 | |
|---|
| 262 | ll1 = MAX( orlgth1, commonAlloc1 ); |
|---|
| 263 | ll2 = MAX( orlgth2, commonAlloc2 ); |
|---|
| 264 | |
|---|
| 265 | fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 ); |
|---|
| 266 | |
|---|
| 267 | commonIP = AllocateIntMtx( ll1+10, ll2+10 ); |
|---|
| 268 | |
|---|
| 269 | fprintf( stderr, "succeeded\n\n" ); |
|---|
| 270 | |
|---|
| 271 | commonAlloc1 = ll1; |
|---|
| 272 | commonAlloc2 = ll2; |
|---|
| 273 | } |
|---|
| 274 | ijp = commonIP; |
|---|
| 275 | |
|---|
| 276 | cpmx_calc( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc ); |
|---|
| 277 | cpmx_calc( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc ); |
|---|
| 278 | |
|---|
| 279 | match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 ); |
|---|
| 280 | match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 ); |
|---|
| 281 | |
|---|
| 282 | if( outgap == 1 ) |
|---|
| 283 | { |
|---|
| 284 | for( i=1; i<lgth1+1; i++ ) |
|---|
| 285 | { |
|---|
| 286 | initverticalw[i] += penalty * 0.5; |
|---|
| 287 | } |
|---|
| 288 | for( j=1; j<lgth2+1; j++ ) |
|---|
| 289 | { |
|---|
| 290 | currentw[j] += penalty * 0.5; |
|---|
| 291 | } |
|---|
| 292 | } |
|---|
| 293 | |
|---|
| 294 | for( j=0; j<lgth2+1; ++j ) |
|---|
| 295 | { |
|---|
| 296 | m[j] = currentw[j-1] + penalty * 0.5; mp[j] = 0; |
|---|
| 297 | } |
|---|
| 298 | |
|---|
| 299 | lastverticalw[0] = currentw[lgth2-1]; |
|---|
| 300 | |
|---|
| 301 | if( outgap ) lasti = lgth1+1; else lasti = lgth1; |
|---|
| 302 | |
|---|
| 303 | for( i=1; i<lasti; i++ ) |
|---|
| 304 | { |
|---|
| 305 | |
|---|
| 306 | floatncpy( previousw, currentw, lgth2+1 ); |
|---|
| 307 | previousw[0] = initverticalw[i-1]; |
|---|
| 308 | |
|---|
| 309 | match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 ); |
|---|
| 310 | currentw[0] = initverticalw[i]; |
|---|
| 311 | |
|---|
| 312 | mi = previousw[0] + penalty * 0.5; mpi = 0; |
|---|
| 313 | for( j=1; j<lgth2+1; j++ ) |
|---|
| 314 | { |
|---|
| 315 | wm = previousw[j-1]; |
|---|
| 316 | ijp[i][j] = 0; |
|---|
| 317 | |
|---|
| 318 | g = penalty * 0.5; |
|---|
| 319 | x = mi + g; |
|---|
| 320 | if( x > wm ) |
|---|
| 321 | { |
|---|
| 322 | wm = x; |
|---|
| 323 | ijp[i][j] = -( j - mpi ); |
|---|
| 324 | } |
|---|
| 325 | g = penalty * 0.5; |
|---|
| 326 | x = previousw[j-1] + g; |
|---|
| 327 | if( mi <= x ) |
|---|
| 328 | { |
|---|
| 329 | mi = x; |
|---|
| 330 | mpi = j-1; |
|---|
| 331 | } |
|---|
| 332 | |
|---|
| 333 | g = penalty * 0.5; |
|---|
| 334 | x = m[j] + g; |
|---|
| 335 | if( x > wm ) |
|---|
| 336 | { |
|---|
| 337 | wm = x; |
|---|
| 338 | ijp[i][j] = +( i - mp[j] ); |
|---|
| 339 | } |
|---|
| 340 | g = penalty * 0.5; |
|---|
| 341 | x = previousw[j-1] + g; |
|---|
| 342 | if( m[j] <= x ) |
|---|
| 343 | { |
|---|
| 344 | m[j] = x; |
|---|
| 345 | mp[j] = i-1; |
|---|
| 346 | } |
|---|
| 347 | currentw[j] += wm; |
|---|
| 348 | } |
|---|
| 349 | lastverticalw[i] = currentw[lgth2-1]; |
|---|
| 350 | } |
|---|
| 351 | /* |
|---|
| 352 | fprintf( stderr, "\n" ); |
|---|
| 353 | for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] ); |
|---|
| 354 | fprintf( stderr, "#####\n" ); |
|---|
| 355 | for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] ); |
|---|
| 356 | fprintf( stderr, "====>" ); |
|---|
| 357 | for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] ); |
|---|
| 358 | for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] ); |
|---|
| 359 | */ |
|---|
| 360 | Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc ); |
|---|
| 361 | |
|---|
| 362 | resultlen = strlen( mseq1[0] ); |
|---|
| 363 | if( alloclen < resultlen || resultlen > N ) |
|---|
| 364 | { |
|---|
| 365 | fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N ); |
|---|
| 366 | ErrorExit( "LENGTH OVER!\n" ); |
|---|
| 367 | } |
|---|
| 368 | |
|---|
| 369 | for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] ); |
|---|
| 370 | for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] ); |
|---|
| 371 | /* |
|---|
| 372 | fprintf( stderr, "\n" ); |
|---|
| 373 | for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] ); |
|---|
| 374 | fprintf( stderr, "#####\n" ); |
|---|
| 375 | for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] ); |
|---|
| 376 | */ |
|---|
| 377 | return( wm ); |
|---|
| 378 | } |
|---|
