6 #define USE_PENALTY_EX 0
10 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
13 char tmpc = s1[0][i1];
17 *match++ = amino_dis[tmpc][*seq2++];
20 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
24 for( j=0; j<lgth2; j++ )
25 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
29 static float Atracking( float *lasthorizontalw, float *lastverticalw,
30 char **seq1, char **seq2,
31 char **mseq1, char **mseq2,
32 float **cpmx1, float **cpmx2,
35 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
38 lgth1 = strlen( seq1[0] );
39 lgth2 = strlen( seq2[0] );
43 for( i=0; i<lgth1; i++ )
45 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
49 for( i=0; i<lgth1+1; i++ )
53 for( j=0; j<lgth2+1; j++ )
55 ijp[0][j] = -( j + 1 );
59 mseq1[0] += lgth1+lgth2;
61 mseq2[0] += lgth1+lgth2;
63 iin = lgth1; jin = lgth2;
64 limk = lgth1+lgth2 + 1;
65 for( k=0; k<limk; k++ )
67 if( ijp[iin][jin] < 0 )
69 ifi = iin-1; jfi = jin+ijp[iin][jin];
71 else if( ijp[iin][jin] > 0 )
73 ifi = iin-ijp[iin][jin]; jfi = jin-1;
77 ifi = iin-1; jfi = jin-1;
82 *--mseq1[0] = seq1[0][ifi+l];
90 *--mseq2[0] = seq2[0][jfi+l];
93 if( iin <= 0 || jin <= 0 ) break;
94 *--mseq1[0] = seq1[0][ifi];
95 *--mseq2[0] = seq2[0][jfi];
103 float G__align11( char **seq1, char **seq2, int alloclen )
104 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
108 int lasti; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
111 float wm; /* int ?????? */
113 float *currentw, *previousw;
114 float fpenalty = (float)penalty;
115 float fpenalty_ex = (float)penalty_ex;
119 float *mjpt, *prept, *curpt;
125 static float *w1, *w2;
127 static float *initverticalw; /* kufuu sureba iranai */
128 static float *lastverticalw; /* kufuu sureba iranai */
132 static float **cpmx1;
133 static float **cpmx2;
134 static int **intwork;
135 static float **floatwork;
136 static int orlgth1 = 0, orlgth2 = 0;
140 mseq1 = AllocateCharMtx( njob, 0 );
141 mseq2 = AllocateCharMtx( njob, 0 );
145 lgth1 = strlen( seq1[0] );
146 lgth2 = strlen( seq2[0] );
148 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
152 if( orlgth1 > 0 && orlgth2 > 0 )
156 FreeFloatVec( match );
157 FreeFloatVec( initverticalw );
158 FreeFloatVec( lastverticalw );
166 FreeFloatMtx( cpmx1 );
167 FreeFloatMtx( cpmx2 );
169 FreeFloatMtx( floatwork );
170 FreeIntMtx( intwork );
173 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
174 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
177 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
180 w1 = AllocateFloatVec( ll2+2 );
181 w2 = AllocateFloatVec( ll2+2 );
182 match = AllocateFloatVec( ll2+2 );
184 initverticalw = AllocateFloatVec( ll1+2 );
185 lastverticalw = AllocateFloatVec( ll1+2 );
187 m = AllocateFloatVec( ll2+2 );
188 mp = AllocateIntVec( ll2+2 );
190 mseq = AllocateCharMtx( njob, ll1+ll2 );
192 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
193 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
195 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
196 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
199 fprintf( stderr, "succeeded\n" );
211 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
215 if( commonAlloc1 && commonAlloc2 )
217 FreeShortMtx( commonIP );
220 ll1 = MAX( orlgth1, commonAlloc1 );
221 ll2 = MAX( orlgth2, commonAlloc2 );
224 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
227 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
230 fprintf( stderr, "succeeded\n\n" );
240 for( i=0; i<lgth1; i++ )
241 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
248 match_calc( initverticalw, seq2, seq1, 0, lgth1 );
251 match_calc( currentw, seq1, seq2, 0, lgth2 );
255 for( i=1; i<lgth1+1; i++ )
257 initverticalw[i] += fpenalty;
259 for( j=1; j<lgth2+1; j++ )
261 currentw[j] += fpenalty;
265 for( j=1; j<lgth2+1; ++j )
267 m[j] = currentw[j-1]; mp[j] = 0;
270 lastverticalw[0] = currentw[lgth2-1];
272 if( outgap ) lasti = lgth1+1; else lasti = lgth1;
275 fprintf( stderr, "currentw = \n" );
276 for( i=0; i<lgth1+1; i++ )
278 fprintf( stderr, "%5.2f ", currentw[i] );
280 fprintf( stderr, "\n" );
281 fprintf( stderr, "initverticalw = \n" );
282 for( i=0; i<lgth2+1; i++ )
284 fprintf( stderr, "%5.2f ", initverticalw[i] );
286 fprintf( stderr, "\n" );
289 for( i=1; i<lasti; i++ )
292 previousw = currentw;
295 previousw[0] = initverticalw[i-1];
297 match_calc( currentw, seq1, seq2, i, lgth2 );
299 fprintf( stderr, "\n" );
300 fprintf( stderr, "i=%d\n", i );
301 fprintf( stderr, "currentw = \n" );
302 for( j=0; j<lgth2; j++ )
304 fprintf( stderr, "%5.2f ", currentw[j] );
306 fprintf( stderr, "\n" );
309 fprintf( stderr, "\n" );
310 fprintf( stderr, "i=%d\n", i );
311 fprintf( stderr, "currentw = \n" );
312 for( j=0; j<lgth2; j++ )
314 fprintf( stderr, "%5.2f ", currentw[j] );
316 fprintf( stderr, "\n" );
318 currentw[0] = initverticalw[i];
320 mi = previousw[0]; mpi = 0;
325 curpt = currentw + 1;
327 for( j=1; j<lgth2+1; j++ )
333 fprintf( stderr, "%5.0f->", wm );
337 fprintf( stderr, "%5.0f?", g );
342 *ijppt = -( j - mpi );
354 g = *mjpt + fpenalty;
356 fprintf( stderr, "%5.0f?", g );
361 *ijppt = +( i - *mpjpt );
374 fprintf( stderr, "%5.0f ", wm );
383 lastverticalw[i] = currentw[lgth2-1];
386 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp );
389 resultlen = strlen( mseq1[0] );
390 if( alloclen < resultlen || resultlen > N )
392 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
393 ErrorExit( "LENGTH OVER!\n" );
397 strcpy( seq1[0], mseq1[0] );
398 strcpy( seq2[0], mseq2[0] );
400 fprintf( stdout, "\n" );
401 fprintf( stdout, ">\n%s\n", mseq1[0] );
402 fprintf( stdout, ">\n%s\n", mseq2[0] );
405 fprintf( stderr, "wm = %f\n", wm );