~ubuntu-branches/debian/sid/mafft/sid : revision 10

1

#include "mltaln.h"

2

#include "dp.h"

3

4

#define MACHIGAI 0

5

#define OUTGAP0TRY 1

6

#define DEBUG 0

7

#define XXXXXXX 0

8

#define USE_PENALTY_EX 0

9

#define FASTMATCHCALC 1

10

11

12

static float **impmtx = NULL;

13

#if 1 // tditeration to naiveQscore_imp de tsukawareru.

14

float imp_match_out_scQ( int i1, int j1 )

15

{

16

// fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );

17

// fprintf( stderr, "val = %f\n", impmtx[i1][j1] );

18

return( impmtx[i1][j1] );

19

}

20

#endif

21

22

static void imp_match_out_veadQ( float *imp, int i1, int lgth2 )

23

{

24

#if FASTMATCHCALC

25

float *pt = impmtx[i1];

26

while( lgth2-- )

27

*imp++ += *pt++;

28

#else

29

int j;

30

float *pt = impmtx[i1];

31

for( j=0; j<lgth2; j++ )

32

*imp++ += pt[j];

33

#endif

34

}

35

static void imp_match_out_vead_tateQ( float *imp, int j1, int lgth1 )

36

{

37

int i;

38

for( i=0; i<lgth1; i++ )

39

*imp++ += impmtx[i][j1];

40

}

41

42

#if 1 // tbfast.c kara yobareru.

43

void imp_match_init_strictQ( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom, int forscore )

44

{

45

int i, j, k1, k2, tmpint, start1, start2, end1, end2;

46

static int impalloclen = 0;

47

float effij;

48

double effijx;

49

char *pt, *pt1, *pt2;

50

static char *nocount1 = NULL;

51

static char *nocount2 = NULL;

52

LocalHom *tmpptr;

53

54

if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )

55

{

56

if( impmtx ) FreeFloatMtx( impmtx );

57

if( nocount1 ) free( nocount1 );

58

if( nocount2 ) free( nocount2 );

59

impalloclen = MAX( lgth1, lgth2 ) + 2;

60

impmtx = AllocateFloatMtx( impalloclen, impalloclen );

61

nocount1 = AllocateCharVec( impalloclen );

62

nocount2 = AllocateCharVec( impalloclen );

63

}

64

65

for( i=0; i<lgth1; i++ )

66

{

67

for( j=0; j<clus1; j++ )

68

if( seq1[j][i] == '-' ) break;

69

if( j != clus1 ) nocount1[i] = 1;

70

else nocount1[i] = 0;

71

}

72

for( i=0; i<lgth2; i++ )

73

{

74

for( j=0; j<clus2; j++ )

75

if( seq2[j][i] == '-' ) break;

76

if( j != clus2 ) nocount2[i] = 1;

77

else nocount2[i] = 0;

78

}

79

80

#if 0

81

fprintf( stderr, "nocount2 =\n" );

82

for( i = 0; i<impalloclen; i++ )

83

{

84

fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );

85

}

86

#endif

87

88

89

90

#if 0

91

fprintf( stderr, "eff1 in _init_strict = \n" );

92

for( i=0; i<clus1; i++ )

93

fprintf( stderr, "eff1[] = %f\n", eff1[i] );

94

for( i=0; i<clus2; i++ )

95

fprintf( stderr, "eff2[] = %f\n", eff2[i] );

96

#endif

97

98

for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )

99

impmtx[i][j] = 0.0;

100

effijx = fastathreshold;

101

for( i=0; i<clus1; i++ )

102

{

103

for( j=0; j<clus2; j++ )

104

{

105

effij = (float)( eff1[i] * eff2[j] * effijx );

106

tmpptr = localhom[i][j];

107

while( tmpptr )

108

{

109

// fprintf( stderr, "start1 = %d\n", tmpptr->start1 );

110

// fprintf( stderr, "end1 = %d\n", tmpptr->end1 );

111

// fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );

112

// fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );

113

pt = seq1[i];

114

tmpint = -1;

115

while( *pt != 0 )

116

{

117

if( *pt++ != '-' ) tmpint++;

118

if( tmpint == tmpptr->start1 ) break;

119

}

120

start1 = pt - seq1[i] - 1;

121

122

if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;

123

else

124

{

125

#if MACHIGAI

126

while( *pt != 0 )

127

{

128

// fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );

129

if( tmpint == tmpptr->end1 ) break;

130

if( *pt++ != '-' ) tmpint++;

131

}

132

end1 = pt - seq1[i] - 0;

133

#else

134

while( *pt != 0 )

135

{

136

// fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );

137

if( *pt++ != '-' ) tmpint++;

138

if( tmpint == tmpptr->end1 ) break;

139

}

140

end1 = pt - seq1[i] - 1;

141

#endif

142

}

143

144

pt = seq2[j];

145

tmpint = -1;

146

while( *pt != 0 )

147

{

148

if( *pt++ != '-' ) tmpint++;

149

if( tmpint == tmpptr->start2 ) break;

150

}

151

start2 = pt - seq2[j] - 1;

152

if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;

153

else

154

{

155

#if MACHIGAI

156

while( *pt != 0 )

157

{

158

if( tmpint == tmpptr->end2 ) break;

159

if( *pt++ != '-' ) tmpint++;

160

}

161

end2 = pt - seq2[j] - 0;

162

#else

163

while( *pt != 0 )

164

{

165

if( *pt++ != '-' ) tmpint++;

166

if( tmpint == tmpptr->end2 ) break;

167

}

168

end2 = pt - seq2[j] - 1;

169

#endif

170

}

171

// fprintf( stderr, "start1 = %d (%c), end1 = %d (%c), start2 = %d (%c), end2 = %d (%c)\n", start1, seq1[i][start1], end1, seq1[i][end1], start2, seq2[j][start2], end2, seq2[j][end2] );

172

// fprintf( stderr, "step 0\n" );

173

if( end1 - start1 != end2 - start2 )

174

{

175

// fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );

176

}

177

178

#if 1

179

k1 = start1; k2 = start2;

180

pt1 = seq1[i] + k1;

181

pt2 = seq2[j] + k2;

182

while( *pt1 && *pt2 )

183

{

184

if( *pt1 != '-' && *pt2 != '-' )

185

{

186

// �Ťߤ��Ťˤ��ʤ��褦��դ��Ʋ��

187

// impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;

188

// impmtx[k1][k2] += tmpptr->importance * effij;

189

impmtx[k1][k2] += tmpptr->fimportance * effij;

190

// fprintf( stderr, "#### impmtx[k1][k2] = %f, tmpptr->fimportance=%f, effij=%f\n", impmtx[k1][k2], tmpptr->fimportance, effij );

191

// fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );

192

// fprintf( stderr, "%d (%c) - %d (%c) - %f\n", k1, *pt1, k2, *pt2, tmpptr->fimportance * effij );

193

k1++; k2++;

194

pt1++; pt2++;

195

}

196

else if( *pt1 != '-' && *pt2 == '-' )

197

{

198

// fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );

199

k2++; pt2++;

200

}

201

else if( *pt1 == '-' && *pt2 != '-' )

202

{

203

// fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );

204

k1++; pt1++;

205

}

206

else if( *pt1 == '-' && *pt2 == '-' )

207

{

208

// fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );

209

k1++; pt1++;

210

k2++; pt2++;

211

}

212

if( k1 > end1 || k2 > end2 ) break;

213

}

214

#else

215

while( k1 <= end1 && k2 <= end2 )

216

{

217

fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 );

218

if( !nocount1[k1] && !nocount2[k2] )

219

{

220

impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j] * fastathreshold;

221

fprintf( stderr, "marked\n" );

222

}

223

else

224

fprintf( stderr, "no count\n" );

225

k1++; k2++;

226

}

227

#endif

228

tmpptr = tmpptr->next;

229

}

230

}

231

}

232

#if 0

233

if( clus1 == 1 && clus2 == 6 )

234

{

235

fprintf( stderr, "\n" );

236

fprintf( stderr, "seq1[0] = %s\n", seq1[0] );

237

fprintf( stderr, "seq2[0] = %s\n", seq2[0] );

238

fprintf( stderr, "impmtx = \n" );

239

for( k2=0; k2<lgth2; k2++ )

240

fprintf( stderr, "%6.3f ", (double)k2 );

241

fprintf( stderr, "\n" );

242

for( k1=0; k1<lgth1; k1++ )

243

{

244

fprintf( stderr, "%d ", k1 );

245

for( k2=0; k2<3; k2++ )

246

fprintf( stderr, "%2.1f ", impmtx[k1][k2] );

247

fprintf( stderr, "\n" );

248

}

249

exit( 1 );

250

}

251

#endif

252

}

253

#endif

254

255

256

static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )

257

{

258

#if FASTMATCHCALC

259

int j, l;

260

float scarr[26];

261

float **cpmxpd = floatwork;

262

int **cpmxpdn = intwork;

263

float *matchpt, *cpmxpdpt, **cpmxpdptpt;

264

int *cpmxpdnpt, **cpmxpdnptpt;

265

if( initialize )

266

{

267

int count = 0;

268

for( j=0; j<lgth2; j++ )

269

{

270

count = 0;

271

for( l=0; l<26; l++ )

272

{

273

if( cpmx2[l][j] )

274

{

275

cpmxpd[j][count] = cpmx2[l][j];

276

cpmxpdn[j][count] = l;

277

count++;

278

}

279

}

280

cpmxpdn[j][count] = -1;

281

}

282

}

283

284

{

285

for( l=0; l<26; l++ )

286

{

287

scarr[l] = 0.0;

288

for( j=0; j<26; j++ )

289

scarr[l] += n_dis[j][l] * cpmx1[j][i1];

290

}

291

matchpt = match;

292

cpmxpdnptpt = cpmxpdn;

293

cpmxpdptpt = cpmxpd;

294

while( lgth2-- )

295

{

296

*matchpt = 0.0;

297

cpmxpdnpt = *cpmxpdnptpt++;

298

cpmxpdpt = *cpmxpdptpt++;

299

while( *cpmxpdnpt>-1 )

300

*matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;

301

matchpt++;

302

}

303

}

304

#else

305

int j, k, l;

306

float scarr[26];

307

float **cpmxpd = floatwork;

308

int **cpmxpdn = intwork;

309

// simple

310

if( initialize )

311

{

312

int count = 0;

313

for( j=0; j<lgth2; j++ )

314

{

315

count = 0;

316

for( l=0; l<26; l++ )

317

{

318

if( cpmx2[l][j] )

319

{

320

cpmxpd[count][j] = cpmx2[l][j];

321

cpmxpdn[count][j] = l;

322

count++;

323

}

324

}

325

cpmxpdn[count][j] = -1;

326

}

327

}

328

for( l=0; l<26; l++ )

329

{

330

scarr[l] = 0.0;

331

for( k=0; k<26; k++ )

332

scarr[l] += n_dis[k][l] * cpmx1[k][i1];

333

}

334

for( j=0; j<lgth2; j++ )

335

{

336

match[j] = 0.0;

337

for( k=0; cpmxpdn[k][j]>-1; k++ )

338

match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];

339

}

340

#endif

341

}

342

343

static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw,

344

char **seq1, char **seq2,

345

char **mseq1, char **mseq2,

346

float **cpmx1, float **cpmx2,

347

short **ijp, int icyc, int jcyc )

348

{

349

int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;

350

float wm;

351

char *gaptable1, *gt1bk;

352

char *gaptable2, *gt2bk;

353

lgth1 = strlen( seq1[0] );

354

lgth2 = strlen( seq2[0] );

355

gt1bk = AllocateCharVec( lgth1+lgth2+1 );

356

gt2bk = AllocateCharVec( lgth1+lgth2+1 );

357

358

#if 0

359

for( i=0; i<lgth1; i++ )

360

{

361

fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );

362

}

363

#endif

364

365

if( outgap == 1 )

366

;

367

else

368

{

369

wm = lastverticalw[0];

370

for( i=0; i<lgth1; i++ )

371

{

372

if( lastverticalw[i] >= wm )

373

{

374

wm = lastverticalw[i];

375

iin = i; jin = lgth2-1;

376

ijp[lgth1][lgth2] = +( lgth1 - i );

377

}

378

}

379

for( j=0; j<lgth2; j++ )

380

{

381

if( lasthorizontalw[j] >= wm )

382

{

383

wm = lasthorizontalw[j];

384

iin = lgth1-1; jin = j;

385

ijp[lgth1][lgth2] = -( lgth2 - j );

386

}

387

}

388

}

389

390

for( i=0; i<lgth1+1; i++ )

391

{

392

ijp[i][0] = i + 1;

393

}

394

for( j=0; j<lgth2+1; j++ )

395

{

396

ijp[0][j] = -( j + 1 );

397

}

398

399

gaptable1 = gt1bk + lgth1+lgth2;

400

*gaptable1 = 0;

401

gaptable2 = gt2bk + lgth1+lgth2;

402

*gaptable2 = 0;

403

404

iin = lgth1; jin = lgth2;

405

*impwmpt = 0.0;

406

for( k=0; k<=lgth1+lgth2; k++ )

407

{

408

if( ijp[iin][jin] < 0 )

409

{

410

ifi = iin-1; jfi = jin+ijp[iin][jin];

411

}

412

else if( ijp[iin][jin] > 0 )

413

{

414

ifi = iin-ijp[iin][jin]; jfi = jin-1;

415

}

416

else

417

{

418

ifi = iin-1; jfi = jin-1;

419

}

420

l = iin - ifi;

421

while( --l )

422

{

423

*--gaptable1 = 'o';

424

*--gaptable2 = '-';

425

k++;

426

}

427

l= jin - jfi;

428

while( --l )

429

{

430

*--gaptable1 = '-';

431

*--gaptable2 = 'o';

432

k++;

433

}

434

if( iin == lgth1 || jin == lgth2 )

435

;

436

else

437

{

438

*impwmpt += imp_match_out_scQ( iin, jin );

439

440

// fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );

441

}

442

if( iin <= 0 || jin <= 0 ) break;

443

*--gaptable1 = 'o';

444

*--gaptable2 = 'o';

445

k++;

446

iin = ifi; jin = jfi;

447

}

448

449

for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );

450

for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );

451

452

free( gt1bk );

453

free( gt2bk );

454

}

455

456

static float Atracking( float *lasthorizontalw, float *lastverticalw,

457

char **seq1, char **seq2,

458

char **mseq1, char **mseq2,

459

float **cpmx1, float **cpmx2,

460

short **ijp, int icyc, int jcyc )

461

{

462

int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;

463

float wm;

464

char *gaptable1, *gt1bk;

465

char *gaptable2, *gt2bk;

466

lgth1 = strlen( seq1[0] );

467

lgth2 = strlen( seq2[0] );

468

469

gt1bk = AllocateCharVec( lgth1+lgth2+1 );

470

gt2bk = AllocateCharVec( lgth1+lgth2+1 );

471

472

#if 0

473

for( i=0; i<lgth1; i++ )

474

{

475

fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );

476

}

477

#endif

478

479

if( outgap == 1 )

480

;

481

else

482

{

483

wm = lastverticalw[0];

484

for( i=0; i<lgth1; i++ )

485

{

486

if( lastverticalw[i] >= wm )

487

{

488

wm = lastverticalw[i];

489

iin = i; jin = lgth2-1;

490

ijp[lgth1][lgth2] = +( lgth1 - i );

491

}

492

}

493

for( j=0; j<lgth2; j++ )

494

{

495

if( lasthorizontalw[j] >= wm )

496

{

497

wm = lasthorizontalw[j];

498

iin = lgth1-1; jin = j;

499

ijp[lgth1][lgth2] = -( lgth2 - j );

500

}

501

}

502

}

503

504

for( i=0; i<lgth1+1; i++ )

505

{

506

ijp[i][0] = i + 1;

507

}

508

for( j=0; j<lgth2+1; j++ )

509

{

510

ijp[0][j] = -( j + 1 );

511

}

512

513

gaptable1 = gt1bk + lgth1+lgth2;

514

*gaptable1 = 0;

515

gaptable2 = gt2bk + lgth1+lgth2;

516

*gaptable2 = 0;

517

518

iin = lgth1; jin = lgth2;

519

for( k=0; k<=lgth1+lgth2; k++ )

520

{

521

if( ijp[iin][jin] < 0 )

522

{

523

ifi = iin-1; jfi = jin+ijp[iin][jin];

524

}

525

else if( ijp[iin][jin] > 0 )

526

{

527

ifi = iin-ijp[iin][jin]; jfi = jin-1;

528

}

529

else

530

{

531

ifi = iin-1; jfi = jin-1;

532

}

533

l = iin - ifi;

534

while( --l )

535

{

536

*--gaptable1 = 'o';

537

*--gaptable2 = '-';

538

k++;

539

}

540

l= jin - jfi;

541

while( --l )

542

{

543

*--gaptable1 = '-';

544

*--gaptable2 = 'o';

545

k++;

546

}

547

if( iin <= 0 || jin <= 0 ) break;

548

*--gaptable1 = 'o';

549

*--gaptable2 = 'o';

550

k++;

551

iin = ifi; jin = jfi;

552

}

553

554

for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );

555

for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );

556

557

free( gt1bk );

558

free( gt2bk );

559

560

return( 0.0 );

561

}

562

563

float Q__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2 )

564

/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */

565

{

566

// int k;

567

register int i, j;

568

int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */

569

int lgth1, lgth2;

570

int resultlen;

571

float wm = 0.0; /* int ?????? */

572

float g;

573

float *currentw, *previousw;

574

// float fpenalty = (float)penalty;

575

#if USE_PENALTY_EX

576

float fpenalty_ex = (float)penalty_ex;

577

#endif

578

#if 1

579

float *wtmp;

580

short *ijppt;

581

float *mjpt, *prept, *curpt;

582

int *mpjpt;

583

#endif

584

static float mi, *m;

585

static short **ijp;

586

static int mpi, *mp;

587

static float *w1, *w2;

588

static float *match;

589

static float *initverticalw; /* kufuu sureba iranai */

590

static float *lastverticalw; /* kufuu sureba iranai */

591

static char **mseq1;

592

static char **mseq2;

593

static char **mseq;

594

static float *digf1;

595

static float *digf2;

596

static float *diaf1;

597

static float *diaf2;

598

static float *gapz1;

599

static float *gapz2;

600

static float *gapf1;

601

static float *gapf2;

602

static float *ogcp1g;

603

static float *ogcp2g;

604

static float *fgcp1g;

605

static float *fgcp2g;

606

static float *og_h_dg_n1_p;

607

static float *og_h_dg_n2_p;

608

static float *fg_h_dg_n1_p;

609

static float *fg_h_dg_n2_p;

610

static float *og_t_fg_h_dg_n1_p;

611

static float *og_t_fg_h_dg_n2_p;

612

static float *fg_t_og_h_dg_n1_p;

613

static float *fg_t_og_h_dg_n2_p;

614

static float *gapz_n1;

615

static float *gapz_n2;

616

static float **cpmx1;

617

static float **cpmx2;

618

static int **intwork;

619

static float **floatwork;

620

static int orlgth1 = 0, orlgth2 = 0;

621

float fpenalty = (float)penalty;

622

float tmppenal;

623

float *fg_t_og_h_dg_n2_p_pt;

624

float *og_t_fg_h_dg_n2_p_pt;

625

float *og_h_dg_n2_p_pt;

626

float *fg_h_dg_n2_p_pt;

627

float *gapz_n2_pt0;

628

float *gapz_n2_pt1;

629

float *fgcp2pt;

630

float *ogcp2pt;

631

float fg_t_og_h_dg_n1_p_va;

632

float og_t_fg_h_dg_n1_p_va;

633

float og_h_dg_n1_p_va;

634

float fg_h_dg_n1_p_va;

635

float gapz_n1_va0;

636

float gapz_n1_va1;

637

float fgcp1va;

638

float ogcp1va;

639

float kyokaipenal;

640

641

642

643

#if 0

644

fprintf( stderr, "#### seq1[0] = %s\n", seq1[0] );

645

fprintf( stderr, "#### seq2[0] = %s\n", seq2[0] );

646

#endif

647

if( orlgth1 == 0 )

648

{

649

mseq1 = AllocateCharMtx( njob, 0 );

650

mseq2 = AllocateCharMtx( njob, 0 );

651

}

652

653

654

lgth1 = strlen( seq1[0] );

655

lgth2 = strlen( seq2[0] );

656

#if 0

657

if( lgth1 == 0 || lgth2 == 0 )

658

{

659

fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );

660

}

661

#endif

662

663

if( lgth1 > orlgth1 || lgth2 > orlgth2 )

664

{

665

int ll1, ll2;

666

667

if( orlgth1 > 0 && orlgth2 > 0 )

668

{

669

FreeFloatVec( w1 );

670

FreeFloatVec( w2 );

671

FreeFloatVec( match );

672

FreeFloatVec( initverticalw );

673

FreeFloatVec( lastverticalw );

674

675

FreeFloatVec( m );

676

FreeIntVec( mp );

677

678

FreeCharMtx( mseq );

679

680

FreeFloatVec( digf1 );

681

FreeFloatVec( digf2 );

682

FreeFloatVec( diaf1 );

683

FreeFloatVec( diaf2 );

684

FreeFloatVec( gapz1 );

685

FreeFloatVec( gapz2 );

686

FreeFloatVec( gapf1 );

687

FreeFloatVec( gapf2 );

688

FreeFloatVec( ogcp1g );

689

FreeFloatVec( ogcp2g );

690

FreeFloatVec( fgcp1g );

691

FreeFloatVec( fgcp2g );

692

FreeFloatVec( og_h_dg_n1_p );

693

FreeFloatVec( og_h_dg_n2_p );

694

FreeFloatVec( fg_h_dg_n1_p );

695

FreeFloatVec( fg_h_dg_n2_p );

696

FreeFloatVec( og_t_fg_h_dg_n1_p );

697

FreeFloatVec( og_t_fg_h_dg_n2_p );

698

FreeFloatVec( fg_t_og_h_dg_n1_p );

699

FreeFloatVec( fg_t_og_h_dg_n2_p );

700

FreeFloatVec( gapz_n1 );

701

FreeFloatVec( gapz_n2 );

702

703

FreeFloatMtx( cpmx1 );

704

FreeFloatMtx( cpmx2 );

705

706

FreeFloatMtx( floatwork );

707

FreeIntMtx( intwork );

708

}

709

710

ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;

711

ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;

712

713

#if DEBUG

714

fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );

715

#endif

716

717

w1 = AllocateFloatVec( ll2+2 );

718

w2 = AllocateFloatVec( ll2+2 );

719

match = AllocateFloatVec( ll2+2 );

720

721

initverticalw = AllocateFloatVec( ll1+2 );

722

lastverticalw = AllocateFloatVec( ll1+2 );

723

724

m = AllocateFloatVec( ll2+2 );

725

mp = AllocateIntVec( ll2+2 );

726

727

mseq = AllocateCharMtx( njob, ll1+ll2 );

728

729

digf1 = AllocateFloatVec( ll1+2 );

730

digf2 = AllocateFloatVec( ll2+2 );

731

diaf1 = AllocateFloatVec( ll1+2 );

732

diaf2 = AllocateFloatVec( ll2+2 );

733

gapz1 = AllocateFloatVec( ll1+2 );

734

gapz2 = AllocateFloatVec( ll2+2 );

735

gapf1 = AllocateFloatVec( ll1+2 );

736

gapf2 = AllocateFloatVec( ll2+2 );

737

ogcp1g = AllocateFloatVec( ll1+2 );

738

ogcp2g = AllocateFloatVec( ll2+2 );

739

fgcp1g = AllocateFloatVec( ll1+2 );

740

fgcp2g = AllocateFloatVec( ll2+2 );

741

og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );

742

og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );

743

fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );

744

fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );

745

og_t_fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );

746

og_t_fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );

747

fg_t_og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 );

748

fg_t_og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 );

749

gapz_n1 = AllocateFloatVec( ll1+2 );

750

gapz_n2 = AllocateFloatVec( ll2+2 );

751

752

cpmx1 = AllocateFloatMtx( 26, ll1+2 );

753

cpmx2 = AllocateFloatMtx( 26, ll2+2 );

754

fprintf( stderr, "ll2 = %d\n", ll2 );

755

756

#if FASTMATCHCALC

757

floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );

758

intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 );

759

#else

760

floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );

761

intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );

762

#endif

763

764

#if DEBUG

765

fprintf( stderr, "succeeded\n" );

766

#endif

767

768

orlgth1 = ll1 - 100;

769

orlgth2 = ll2 - 100;

770

}

771

772

773

for( i=0; i<icyc; i++ )

774

{

775

mseq1[i] = mseq[i];

776

seq1[i][lgth1] = 0;

777

}

778

for( j=0; j<jcyc; j++ )

779

{

780

mseq2[j] = mseq[icyc+j];

781

seq2[j][lgth2] = 0;

782

}

783

784

785

if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )

786

{

787

int ll1, ll2;

788

789

if( commonAlloc1 && commonAlloc2 )

790

{

791

FreeShortMtx( commonIP );

792

}

793

794

ll1 = MAX( orlgth1, commonAlloc1 );

795

ll2 = MAX( orlgth2, commonAlloc2 );

796

797

#if DEBUG

798

fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );

799

#endif

800

801

commonIP = AllocateShortMtx( ll1+10, ll2+10 );

802

803

#if DEBUG

804

fprintf( stderr, "succeeded\n\n" );

805

#endif

806

807

commonAlloc1 = ll1;

808

commonAlloc2 = ll2;

809

}

810

ijp = commonIP;

811

812

#if 0

813

{

814

float t = 0.0;

815

for( i=0; i<icyc; i++ )

816

t += eff1[i];

817

fprintf( stderr, "## totaleff = %f\n", t );

818

}

819

#endif

820

821

cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );

822

cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );

823

824

if( sgap1 )

825

{

826

new_OpeningGapCount_zure( ogcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );

827

new_OpeningGapCount_zure( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );

828

new_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 );

829

new_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 );

830

getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );

831

getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );

832

getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 );

833

getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 );

834

getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );

835

getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );

836

getgapfreq_zure_part( gapz1, icyc, seq1, eff1, lgth1, sgap1 );

837

getgapfreq_zure_part( gapz2, jcyc, seq2, eff2, lgth2, sgap1 );

838

}

839

else

840

{

841

st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 );

842

st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 );

843

st_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1 );

844

st_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2 );

845

getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 );

846

getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 );

847

getdiaminofreq_x( diaf1, icyc, seq1, eff1, lgth1 );

848

getdiaminofreq_x( diaf2, jcyc, seq2, eff2, lgth2 );

849

getgapfreq( gapf1, icyc, seq1, eff1, lgth1 );

850

getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 );

851

getgapfreq_zure( gapz1, icyc, seq1, eff1, lgth1 );

852

getgapfreq_zure( gapz2, jcyc, seq2, eff2, lgth2 );

853

}

854

855

#if 1 // wakaran

856

lastj = lgth2+2;

857

// fprintf( stderr, "lgth2 = %d\n", lgth2 );

858

for( i=0; i<lastj; i++ )

859

{

860

og_h_dg_n2_p[i] = ( 1.0-ogcp2g[i]-digf2[i] ) * fpenalty * 0.5;

861

fg_h_dg_n2_p[i] = ( 1.0-fgcp2g[i]-digf2[i] ) * fpenalty * 0.5;

862

og_t_fg_h_dg_n2_p[i] = (1.0-ogcp2g[i]+fgcp2g[i]-digf2[i]) * 0.5 * fpenalty;

863

fg_t_og_h_dg_n2_p[i] = (1.0-fgcp2g[i]+ogcp2g[i]-digf2[i]) * 0.5 * fpenalty;

864

gapz_n2[i] = (1.0-gapz2[i]);

865

}

866

lastj = lgth1+2;

867

for( i=0; i<lastj; i++ )

868

{

869

og_h_dg_n1_p[i] = ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;

870

fg_h_dg_n1_p[i] = ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;

871

og_t_fg_h_dg_n1_p[i] = (1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) * 0.5 * fpenalty;

872

fg_t_og_h_dg_n1_p[i] = (1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) * 0.5 * fpenalty;

873

gapz_n1[i] = (1.0-gapz1[i]);

874

}

875

#endif

876

877

878

#if 0

879

for( i=0; i<lgth1; i++ )

880

fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );

881

#endif

882

883

currentw = w1;

884

previousw = w2;

885

886

match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );

887

if( localhom )

888

imp_match_out_vead_tateQ( initverticalw, 0, lgth1 ); // 060306

889

890

match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );

891

if( localhom )

892

imp_match_out_veadQ( currentw, 0, lgth2 ); // 060306

893

894

895

#if 0 // -> tbfast.c

896

if( localhom )

897

imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );

898

899

#endif

900

901

902

903

kyokaipenal = 0.0;

904

if( outgap == 1 )

905

{

906

g = 0.0;

907

908

// g += ogcp1g[0] * og_h_dg_n2_p[0];

909

g += ogcp1g[0] * ( 1.0-ogcp2g[0]-digf2[0] ) * fpenalty * 0.5;

910

// if( g ) fprintf( stderr, "init-match penal1=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );

911

912

// g += ogcp2g[0] * og_h_dg_n1_p[0];

913

g += ogcp2g[0] * ( 1.0-ogcp1g[0]-digf1[0] ) * fpenalty * 0.5;

914

// if( g ) fprintf( stderr, "init-match penal2=%f, %c-%c\n", g, seq1[0][0], seq2[0][0] );

915

916

// g += fgcp1g[0] * fg_h_dg_n2_p[0];

917

g += fgcp1g[0] * ( 1.0-fgcp2g[0]-digf2[0] ) * fpenalty * 0.5;

918

// if( g ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );

919

920

// g += fgcp2g[0] * fg_h_dg_n1_p[0];

921

g += fgcp2g[0] * ( 1.0-fgcp1g[0]-digf1[0] ) * fpenalty * 0.5;

922

// if( g ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );

923

924

kyokaipenal = g;

925

initverticalw[0] += g;

926

currentw[0] += g;

927

928

for( i=1; i<lgth1+1; i++ )

929

{

930

// tmppenal = gapz_n2[0]*og_t_fg_h_dg_n1_p[0];

931

tmppenal = ( (1.0-gapz2[0])*(1.0-ogcp1g[0]+fgcp1g[0]-digf1[0]) ) * 0.5 * fpenalty; // mada

932

initverticalw[i] += tmppenal;

933

934

// tmppenal = gapz_n2[1]*fg_t_og_h_dg_n1_p[i];

935

tmppenal = ( (1.0-gapz2[1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada

936

initverticalw[i] += tmppenal;

937

938

}

939

for( j=1; j<lgth2+1; j++ )

940

{

941

// tmppenal = gapz_n1[0]*og_t_fg_h_dg_n2_p[0];

942

tmppenal = ( (1.0-gapz1[0])*(1.0-ogcp2g[0]+fgcp2g[0]-digf2[0]) ) * 0.5 * fpenalty; // mada

943

currentw[j] += tmppenal;

944

945

// tmppenal = gapz_n1[1]*fg_t_og_h_dg_n2_p[j];

946

tmppenal = ( (1.0-gapz1[1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada

947

currentw[j] += tmppenal;

948

}

949

}

950

#if OUTGAP0TRY

951

else

952

{

953

for( j=1; j<lgth2+1; j++ )

954

currentw[j] -= offset * j / 2.0;

955

for( i=1; i<lgth1+1; i++ )

956

initverticalw[i] -= offset * i / 2.0;

957

}

958

#endif

959

960

m[0] = 0.0; // iranai

961

for( j=1; j<lgth2+1; ++j )

962

{

963

mp[j] = 0;

964

m[j] = currentw[j-1] + 10000 * fpenalty; //iinoka?

965

}

966

if( lgth2 == 0 )

967

lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari

968

else

969

lastverticalw[0] = currentw[lgth2-1];

970

971

if( outgap ) lasti = lgth1+1; else lasti = lgth1;

972

973

#if XXXXXXX

974

fprintf( stderr, "currentw = \n" );

975

for( i=0; i<lgth1+1; i++ )

976

{

977

fprintf( stderr, "%5.2f ", currentw[i] );

978

}

979

fprintf( stderr, "\n" );

980

fprintf( stderr, "initverticalw = \n" );

981

for( i=0; i<lgth2+1; i++ )

982

{

983

fprintf( stderr, "%5.2f ", initverticalw[i] );

984

}

985

fprintf( stderr, "\n" );

986

fprintf( stderr, "fcgp\n" );

987

for( i=0; i<lgth1; i++ )

988

fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );

989

for( i=0; i<lgth2; i++ )

990

fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );

991

#endif

992

993

for( i=1; i<lasti; i++ )

994

{

995

wtmp = previousw;

996

previousw = currentw;

997

currentw = wtmp;

998

999

previousw[0] = initverticalw[i-1];

1000

1001

match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );

1002

#if XXXXXXX

1003

fprintf( stderr, "\n" );

1004

fprintf( stderr, "i=%d\n", i );

1005

fprintf( stderr, "currentw = \n" );

1006

for( j=0; j<lgth2; j++ )

1007

{

1008

fprintf( stderr, "%5.2f ", currentw[j] );

1009

}

1010

fprintf( stderr, "\n" );

1011

#endif

1012

if( localhom )

1013

{

1014

// fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );

1015

#if 0

1016

imp_match_out_veadQ( currentw, i, lgth2 );

1017

#else

1018

imp_match_out_veadQ( currentw, i, lgth2 );

1019

#endif

1020

}

1021

#if XXXXXXX

1022

fprintf( stderr, "\n" );

1023

fprintf( stderr, "i=%d\n", i );

1024

fprintf( stderr, "currentw = \n" );

1025

for( j=0; j<lgth2; j++ )

1026

{

1027

fprintf( stderr, "%5.2f ", currentw[j] );

1028

}

1029

fprintf( stderr, "\n" );

1030

#endif

1031

currentw[0] = initverticalw[i];

1032

1033

mpi = 0;

1034

mi = previousw[0] + 10000 * fpenalty;

1035

1036

ijppt = ijp[i] + 1;

1037

mjpt = m + 1;

1038

prept = previousw;

1039

curpt = currentw + 1;

1040

mpjpt = mp + 1;

1041

fg_t_og_h_dg_n2_p_pt = fg_t_og_h_dg_n2_p + 1;

1042

og_t_fg_h_dg_n2_p_pt = og_t_fg_h_dg_n2_p + 1;

1043

og_h_dg_n2_p_pt = og_h_dg_n2_p + 1;

1044

fg_h_dg_n2_p_pt = fg_h_dg_n2_p + 1;

1045

gapz_n2_pt0 = gapz_n2 + 1;

1046

gapz_n2_pt1 = gapz_n2 + 2;

1047

fgcp2pt = fgcp2g + 1;

1048

ogcp2pt = ogcp2g + 1;

1049

1050

fg_t_og_h_dg_n1_p_va = fg_t_og_h_dg_n1_p[i];

1051

og_t_fg_h_dg_n1_p_va = og_t_fg_h_dg_n1_p[i];

1052

og_h_dg_n1_p_va = og_h_dg_n1_p[i];

1053

fg_h_dg_n1_p_va = fg_h_dg_n1_p[i];

1054

gapz_n1_va0 = gapz_n1[i];

1055

gapz_n1_va1 = gapz_n1[i+1];

1056

fgcp1va = fgcp1g[i];

1057

ogcp1va = ogcp1g[i];

1058

1059

lastj = lgth2+1;

1060

for( j=1; j<lastj; j++ )

1061

{

1062

wm = *prept;

1063

1064

// g = ogcp1va * *og_h_dg_n2_p_pt;

1065

// g = ogcp1g[i] * og_h_dg_n2_p[j];

1066

g = ogcp1g[i] * ( 1.0-ogcp2g[j]-digf2[j] ) * fpenalty * 0.5;

1067

// if( g && i==j ) fprintf( stderr, "match penal1=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );

1068

wm += g;

1069

1070

// g = *ogcp2pt * og_h_dg_n1_p_va;

1071

// g = ogcp2g[j] * og_h_dg_n1_p[i];

1072

g = ogcp2g[j] * ( 1.0-ogcp1g[i]-digf1[i] ) * fpenalty * 0.5;

1073

// if( g && i==j ) fprintf( stderr, "match penal2=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );

1074

wm += g;

1075

1076

// g = fgcp1va * *fg_h_dg_n2_p_pt;

1077

// g = fgcp1g[i] * fg_h_dg_n2_p[j];

1078

g = fgcp1g[i] * ( 1.0-fgcp2g[j]-digf2[j] ) * fpenalty * 0.5;

1079

// if( g && i==j ) fprintf( stderr, "match penal3=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );

1080

wm += g;

1081

1082

// g = *fgcp2pt * fg_h_dg_n1_p_va;

1083

// g = fgcp2g[j] * fg_h_dg_n1_p[i];

1084

g = fgcp2g[j] * ( 1.0-fgcp1g[i]-digf1[i] ) * fpenalty * 0.5;

1085

// if( g && i==j ) fprintf( stderr, "match penal4=%f, %c-%c\n", g, seq1[0][i], seq2[0][j] );

1086

wm += g;

1087

1088

*ijppt = 0;

1089

1090

#if 0

1091

fprintf( stderr, "%5.0f->", wm );

1092

#endif

1093

#if 0

1094

fprintf( stderr, "%5.0f?", g );

1095

#endif

1096

// tmppenal = gapz_n1_va1 * *fg_t_og_h_dg_n2_p_pt;

1097

// tmppenal = gapz_n1[i+1] * fg_t_og_h_dg_n2_p[j];

1098

tmppenal = ( (1.0-gapz1[i+1])*(1.0-fgcp2g[j]+ogcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada

1099

if( (g=mi+tmppenal) > wm )

1100

{

1101

// fprintf( stderr, "jump i start=%f (j=%d, fgcp2g[j]=%f, digf2[j]=%f, diaf2[j]=%f), %c-%c\n", g-mi, j, fgcp2g[j], digf2[j], diaf2[j], seq1[0][i], seq2[0][j] );

1102

wm = g;

1103

*ijppt = -( j - mpi );

1104

}

1105

// tmppenal = gapz_n1_va0 * *og_t_fg_h_dg_n2_p_pt;

1106

// tmppenal = gapz_n1[i] * og_t_fg_h_dg_n2_p[j];

1107

tmppenal = ( (1.0-gapz1[i])*(1.0-ogcp2g[j]+fgcp2g[j]-digf2[j]) ) * 0.5 * fpenalty; // mada

1108

if( (g=*prept+tmppenal) >= mi )

1109

{

1110

// fprintf( stderr, "jump i end=%f, %c-%c\n", g-*prept, seq1[0][i-1], seq2[0][j-1] );

1111

mi = g;

1112

mpi = j-1;

1113

}

1114

1115

#if 0

1116

fprintf( stderr, "%5.0f?", g );

1117

#endif

1118

// tmppenal = *gapz_n2_pt1 * fg_t_og_h_dg_n1_p_va;

1119

// tmppenal = gapz_n2[j+1] * fg_t_og_h_dg_n1_p[i];

1120

tmppenal = ( (1.0-gapz2[j+1])*(1.0-fgcp1g[i]+ogcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada

1121

if( (g=*mjpt+tmppenal) > wm )

1122

{

1123

wm = g;

1124

*ijppt = +( i - *mpjpt );

1125

}

1126

// tmppenal = *gapz_n2_pt0 * og_t_fg_h_dg_n1_p_va;

1127

// tmppenal = gapz_n2[j] * og_t_fg_h_dg_n1_p[i];

1128

tmppenal = ( (1.0-gapz2[j])*(1.0-ogcp1g[i]+fgcp1g[i]-digf1[i]) ) * 0.5 * fpenalty; // mada

1129

if( (g=*prept+tmppenal) >= *mjpt )

1130

{

1131

*mjpt = g;

1132

*mpjpt = i-1;

1133

}

1134

#if 0

1135

fprintf( stderr, "%5.0f ", wm );

1136

#endif

1137

*curpt++ += wm;

1138

ijppt++;

1139

mjpt++;

1140

prept++;

1141

mpjpt++;

1142

fg_t_og_h_dg_n2_p_pt++;

1143

og_t_fg_h_dg_n2_p_pt++;

1144

og_h_dg_n2_p_pt++;

1145

fg_h_dg_n2_p_pt++;

1146

gapz_n2_pt0++;

1147

gapz_n2_pt1++;

1148

fgcp2pt++;

1149

ogcp2pt++;

1150

}

1151

lastverticalw[i] = currentw[lgth2-1];

1152

}

1153

1154

// fprintf( stderr, "wm = %f\n", wm );

1155

1156

#if OUTGAP0TRY

1157

if( !outgap )

1158

{

1159

for( j=1; j<lgth2+1; j++ )

1160

currentw[j] -= offset * ( lgth2 - j ) / 2.0;

1161

for( i=1; i<lgth1+1; i++ )

1162

lastverticalw[i] -= offset * ( lgth1 - i / 2.0);

1163

}

1164

#endif

1165

1166

/*

1167

fprintf( stderr, "\n" );

1168

for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );

1169

fprintf( stderr, "#####\n" );

1170

for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );

1171

fprintf( stderr, "====>" );

1172

for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );

1173

for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );

1174

*/

1175

if( localhom )

1176

{

1177

Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );

1178

}

1179

else

1180

Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );

1181

1182

// fprintf( stderr, "### impmatch = %f\n", *impmatch );

1183

1184

resultlen = strlen( mseq1[0] );

1185

if( alloclen < resultlen || resultlen > N )

1186

{

1187

fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );

1188

ErrorExit( "LENGTH OVER!\n" );

1189

}

1190

1191

1192

for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );

1193

for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );

1194

/*

1195

fprintf( stderr, "\n" );

1196

for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );

1197

fprintf( stderr, "#####\n" );

1198

for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );

1199

*/

1200

1201

fprintf( stderr, "wm = %f\n", wm );

1202

1203

1204

if( sgap1 )

1205

return( wm );

1206

else

1207

return( wm );

1208

}

1209

1210