1 /* Copyright 2018-2019 by Michiel de Hoon. All rights reserved.
2 * This file is part of the Biopython distribution and governed by your
3 * choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
4 * Please see the LICENSE file that should have been included as part of this
5 * package.
6 */
7
8
9
10 #define PY_SSIZE_T_CLEAN
11 #include "Python.h"
12 #include "float.h"
13
14
15 #define HORIZONTAL 0x1
16 #define VERTICAL 0x2
17 #define DIAGONAL 0x4
18 #define STARTPOINT 0x8
19 #define ENDPOINT 0x10
20 #define M_MATRIX 0x1
21 #define Ix_MATRIX 0x2
22 #define Iy_MATRIX 0x4
23 #define DONE 0x3
24 #define NONE 0x7
25
26 #define OVERFLOW_ERROR -1
27 #define MEMORY_ERROR -2
28
29 #define MISSING_LETTER -1
30
31 #define SAFE_ADD(t, s) \
32 { if (s != OVERFLOW_ERROR) { \
33 term = t; \
34 if (term > PY_SSIZE_T_MAX - s) s = OVERFLOW_ERROR; \
35 else s += term; \
36 } \
37 }
38
39
40 typedef enum {NeedlemanWunschSmithWaterman,
41 Gotoh,
42 WatermanSmithBeyer,
43 Unknown} Algorithm;
44
45 typedef enum {Global, Local} Mode;
46
47 typedef struct {
48 unsigned char trace : 5;
49 unsigned char path : 3;
50 } Trace;
51
52 typedef struct {
53 unsigned char Ix : 4;
54 unsigned char Iy : 4;
55 } TraceGapsGotoh;
56
57 typedef struct {
58 int* MIx;
59 int* IyIx;
60 int* MIy;
61 int* IxIy;
62 } TraceGapsWatermanSmithBeyer;
63
64 typedef struct {
65 PyObject_HEAD
66 Trace** M;
67 union { TraceGapsGotoh** gotoh;
68 TraceGapsWatermanSmithBeyer** waterman_smith_beyer; } gaps;
69 int nA;
70 int nB;
71 int iA;
72 int iB;
73 Mode mode;
74 Algorithm algorithm;
75 Py_ssize_t length;
76 unsigned char strand;
77 } PathGenerator;
78
79 static PyObject*
PathGenerator_create_path(PathGenerator * self,int i,int j)80 PathGenerator_create_path(PathGenerator* self, int i, int j) {
81 PyObject* tuple;
82 PyObject* row;
83 PyObject* value;
84 int path;
85 const int ii = i;
86 const int jj = j;
87 int n = 1;
88 int direction = 0;
89 Trace** M = self->M;
90 const unsigned char strand = self->strand;
91
92 while (1) {
93 path = M[i][j].path;
94 if (!path) break;
95 if (path != direction) {
96 n++;
97 direction = path;
98 }
99 switch (path) {
100 case HORIZONTAL: j++; break;
101 case VERTICAL: i++; break;
102 case DIAGONAL: i++; j++; break;
103 }
104 }
105
106 i = ii;
107 j = jj;
108 direction = 0;
109 tuple = PyTuple_New(n);
110 if (!tuple) return NULL;
111
112 n = 0;
113 switch (strand) {
114 case '+':
115 while (1) {
116 path = M[i][j].path;
117 if (path != direction) {
118 row = PyTuple_New(2);
119 if (!row) break;
120 value = PyLong_FromLong(i);
121 if (!value) {
122 Py_DECREF(row); /* all references were stolen */
123 break;
124 }
125 PyTuple_SET_ITEM(row, 0, value);
126 value = PyLong_FromLong(j);
127 if (!value) {
128 Py_DECREF(row); /* all references were stolen */
129 break;
130 }
131 PyTuple_SET_ITEM(row, 1, value);
132 PyTuple_SET_ITEM(tuple, n, row);
133 n++;
134 direction = path;
135 }
136 switch (path) {
137 case HORIZONTAL: j++; break;
138 case VERTICAL: i++; break;
139 case DIAGONAL: i++; j++; break;
140 default: return tuple;
141 }
142 }
143 break;
144 case '-': {
145 const int nB = self->nB;
146 while (1) {
147 path = M[i][j].path;
148 if (path != direction) {
149 row = PyTuple_New(2);
150 if (!row) break;
151 value = PyLong_FromLong(i);
152 if (!value) {
153 Py_DECREF(row); /* all references were stolen */
154 break;
155 }
156 PyTuple_SET_ITEM(row, 0, value);
157 value = PyLong_FromLong(nB-j);
158 if (!value) {
159 Py_DECREF(row); /* all references were stolen */
160 break;
161 }
162 PyTuple_SET_ITEM(row, 1, value);
163 PyTuple_SET_ITEM(tuple, n, row);
164 n++;
165 direction = path;
166 }
167 switch (path) {
168 case HORIZONTAL: j++; break;
169 case VERTICAL: i++; break;
170 case DIAGONAL: i++; j++; break;
171 default: return tuple;
172 }
173 }
174 break;
175 }
176 }
177 Py_DECREF(tuple); /* all references were stolen */
178 return PyErr_NoMemory();
179 }
180
181 static Py_ssize_t
PathGenerator_needlemanwunsch_length(PathGenerator * self)182 PathGenerator_needlemanwunsch_length(PathGenerator* self)
183 {
184 int i;
185 int j;
186 int trace;
187 const int nA = self->nA;
188 const int nB = self->nB;
189 Trace** M = self->M;
190 Py_ssize_t term;
191 Py_ssize_t count = MEMORY_ERROR;
192 Py_ssize_t temp;
193 Py_ssize_t* counts;
194 counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
195 if (!counts) goto exit;
196 counts[0] = 1;
197 for (j = 1; j <= nB; j++) {
198 trace = M[0][j].trace;
199 count = 0;
200 if (trace & HORIZONTAL) SAFE_ADD(counts[j-1], count);
201 counts[j] = count;
202 }
203 for (i = 1; i <= nA; i++) {
204 trace = M[i][0].trace;
205 count = 0;
206 if (trace & VERTICAL) SAFE_ADD(counts[0], count);
207 temp = counts[0];
208 counts[0] = count;
209 for (j = 1; j <= nB; j++) {
210 trace = M[i][j].trace;
211 count = 0;
212 if (trace & HORIZONTAL) SAFE_ADD(counts[j-1], count);
213 if (trace & VERTICAL) SAFE_ADD(counts[j], count);
214 if (trace & DIAGONAL) SAFE_ADD(temp, count);
215 temp = counts[j];
216 counts[j] = count;
217 }
218 }
219 PyMem_Free(counts);
220 exit:
221 return count;
222 }
223
224 static Py_ssize_t
PathGenerator_smithwaterman_length(PathGenerator * self)225 PathGenerator_smithwaterman_length(PathGenerator* self)
226 {
227 int i;
228 int j;
229 int trace;
230 const int nA = self->nA;
231 const int nB = self->nB;
232 Trace** M = self->M;
233 Py_ssize_t term;
234 Py_ssize_t count = MEMORY_ERROR;
235 Py_ssize_t total = 0;
236 Py_ssize_t temp;
237 Py_ssize_t* counts;
238 counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
239 if (!counts) goto exit;
240 counts[0] = 1;
241 for (j = 1; j <= nB; j++) counts[j] = 1;
242 for (i = 1; i <= nA; i++) {
243 temp = counts[0];
244 counts[0] = 1;
245 for (j = 1; j <= nB; j++) {
246 trace = M[i][j].trace;
247 count = 0;
248 if (trace & DIAGONAL) SAFE_ADD(temp, count);
249 if (M[i][j].trace & ENDPOINT) SAFE_ADD(count, total);
250 if (trace & HORIZONTAL) SAFE_ADD(counts[j-1], count);
251 if (trace & VERTICAL) SAFE_ADD(counts[j], count);
252 temp = counts[j];
253 if (count == 0 && (trace & STARTPOINT)) count = 1;
254 counts[j] = count;
255 }
256 }
257 count = total;
258 PyMem_Free(counts);
259 exit:
260 return count;
261 }
262
263 static Py_ssize_t
PathGenerator_gotoh_global_length(PathGenerator * self)264 PathGenerator_gotoh_global_length(PathGenerator* self)
265 {
266 int i;
267 int j;
268 int trace;
269 const int nA = self->nA;
270 const int nB = self->nB;
271 Trace** M = self->M;
272 TraceGapsGotoh** gaps = self->gaps.gotoh;
273 Py_ssize_t count = MEMORY_ERROR;
274 Py_ssize_t term;
275 Py_ssize_t M_temp;
276 Py_ssize_t Ix_temp;
277 Py_ssize_t Iy_temp;
278 Py_ssize_t* M_counts = NULL;
279 Py_ssize_t* Ix_counts = NULL;
280 Py_ssize_t* Iy_counts = NULL;
281 M_counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
282 if (!M_counts) goto exit;
283 Ix_counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
284 if (!Ix_counts) goto exit;
285 Iy_counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
286 if (!Iy_counts) goto exit;
287 M_counts[0] = 1;
288 Ix_counts[0] = 0;
289 Iy_counts[0] = 0;
290 for (j = 1; j <= nB; j++) {
291 M_counts[j] = 0;
292 Ix_counts[j] = 0;
293 Iy_counts[j] = 1;
294 }
295 for (i = 1; i <= nA; i++) {
296 M_temp = M_counts[0];
297 M_counts[0] = 0;
298 Ix_temp = Ix_counts[0];
299 Ix_counts[0] = 1;
300 Iy_temp = Iy_counts[0];
301 Iy_counts[0] = 0;
302 for (j = 1; j <= nB; j++) {
303 count = 0;
304 trace = M[i][j].trace;
305 if (trace & M_MATRIX) SAFE_ADD(M_temp, count);
306 if (trace & Ix_MATRIX) SAFE_ADD(Ix_temp, count);
307 if (trace & Iy_MATRIX) SAFE_ADD(Iy_temp, count);
308 M_temp = M_counts[j];
309 M_counts[j] = count;
310 count = 0;
311 trace = gaps[i][j].Ix;
312 if (trace & M_MATRIX) SAFE_ADD(M_temp, count);
313 if (trace & Ix_MATRIX) SAFE_ADD(Ix_counts[j], count);
314 if (trace & Iy_MATRIX) SAFE_ADD(Iy_counts[j], count);
315 Ix_temp = Ix_counts[j];
316 Ix_counts[j] = count;
317 count = 0;
318 trace = gaps[i][j].Iy;
319 if (trace & M_MATRIX) SAFE_ADD(M_counts[j-1], count);
320 if (trace & Ix_MATRIX) SAFE_ADD(Ix_counts[j-1], count);
321 if (trace & Iy_MATRIX) SAFE_ADD(Iy_counts[j-1], count);
322 Iy_temp = Iy_counts[j];
323 Iy_counts[j] = count;
324 }
325 }
326 count = 0;
327 if (M[nA][nB].trace) SAFE_ADD(M_counts[nB], count);
328 if (gaps[nA][nB].Ix) SAFE_ADD(Ix_counts[nB], count);
329 if (gaps[nA][nB].Iy) SAFE_ADD(Iy_counts[nB], count);
330 exit:
331 if (M_counts) PyMem_Free(M_counts);
332 if (Ix_counts) PyMem_Free(Ix_counts);
333 if (Iy_counts) PyMem_Free(Iy_counts);
334 return count;
335 }
336
337 static Py_ssize_t
PathGenerator_gotoh_local_length(PathGenerator * self)338 PathGenerator_gotoh_local_length(PathGenerator* self)
339 {
340 int i;
341 int j;
342 int trace;
343 const int nA = self->nA;
344 const int nB = self->nB;
345 Trace** M = self->M;
346 TraceGapsGotoh** gaps = self->gaps.gotoh;
347 Py_ssize_t term;
348 Py_ssize_t count = MEMORY_ERROR;
349 Py_ssize_t total = 0;
350 Py_ssize_t M_temp;
351 Py_ssize_t Ix_temp;
352 Py_ssize_t Iy_temp;
353 Py_ssize_t* M_counts = NULL;
354 Py_ssize_t* Ix_counts = NULL;
355 Py_ssize_t* Iy_counts = NULL;
356 M_counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
357 if (!M_counts) goto exit;
358 Ix_counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
359 if (!Ix_counts) goto exit;
360 Iy_counts = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
361 if (!Iy_counts) goto exit;
362 M_counts[0] = 1;
363 Ix_counts[0] = 0;
364 Iy_counts[0] = 0;
365 for (j = 1; j <= nB; j++) {
366 M_counts[j] = 1;
367 Ix_counts[j] = 0;
368 Iy_counts[j] = 0;
369 }
370 for (i = 1; i <= nA; i++) {
371 M_temp = M_counts[0];
372 M_counts[0] = 1;
373 Ix_temp = Ix_counts[0];
374 Ix_counts[0] = 0;
375 Iy_temp = Iy_counts[0];
376 Iy_counts[0] = 0;
377 for (j = 1; j <= nB; j++) {
378 count = 0;
379 trace = M[i][j].trace;
380 if (trace & M_MATRIX) SAFE_ADD(M_temp, count);
381 if (trace & Ix_MATRIX) SAFE_ADD(Ix_temp, count);
382 if (trace & Iy_MATRIX) SAFE_ADD(Iy_temp, count);
383 if (count == 0 && (trace & STARTPOINT)) count = 1;
384 M_temp = M_counts[j];
385 M_counts[j] = count;
386 if (M[i][j].trace & ENDPOINT) SAFE_ADD(count, total);
387 count = 0;
388 trace = gaps[i][j].Ix;
389 if (trace & M_MATRIX) SAFE_ADD(M_temp, count);
390 if (trace & Ix_MATRIX) SAFE_ADD(Ix_counts[j], count);
391 if (trace & Iy_MATRIX) SAFE_ADD(Iy_counts[j], count);
392 Ix_temp = Ix_counts[j];
393 Ix_counts[j] = count;
394 count = 0;
395 trace = gaps[i][j].Iy;
396 if (trace & M_MATRIX) SAFE_ADD(M_counts[j-1], count);
397 if (trace & Ix_MATRIX) SAFE_ADD(Ix_counts[j-1], count);
398 if (trace & Iy_MATRIX) SAFE_ADD(Iy_counts[j-1], count);
399 Iy_temp = Iy_counts[j];
400 Iy_counts[j] = count;
401 }
402 }
403 count = total;
404 exit:
405 if (M_counts) PyMem_Free(M_counts);
406 if (Ix_counts) PyMem_Free(Ix_counts);
407 if (Iy_counts) PyMem_Free(Iy_counts);
408 return count;
409 }
410
411 static Py_ssize_t
PathGenerator_waterman_smith_beyer_global_length(PathGenerator * self)412 PathGenerator_waterman_smith_beyer_global_length(PathGenerator* self)
413 {
414 int i;
415 int j;
416 int trace;
417 int* p;
418 int gap;
419 const int nA = self->nA;
420 const int nB = self->nB;
421 Trace** M = self->M;
422 TraceGapsWatermanSmithBeyer** gaps = self->gaps.waterman_smith_beyer;
423 Py_ssize_t count = MEMORY_ERROR;
424 Py_ssize_t term;
425 Py_ssize_t** M_count = NULL;
426 Py_ssize_t** Ix_count = NULL;
427 Py_ssize_t** Iy_count = NULL;
428 M_count = PyMem_Malloc((nA+1)*sizeof(Py_ssize_t*));
429 if (!M_count) goto exit;
430 Ix_count = PyMem_Malloc((nA+1)*sizeof(Py_ssize_t*));
431 if (!Ix_count) goto exit;
432 Iy_count = PyMem_Malloc((nA+1)*sizeof(Py_ssize_t*));
433 if (!Iy_count) goto exit;
434 for (i = 0; i <= nA; i++) {
435 M_count[i] = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
436 if (!M_count[i]) goto exit;
437 Ix_count[i] = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
438 if (!Ix_count[i]) goto exit;
439 Iy_count[i] = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
440 if (!Iy_count[i]) goto exit;
441 }
442 for (i = 0; i <= nA; i++) {
443 for (j = 0; j <= nB; j++) {
444 count = 0;
445 trace = M[i][j].trace;
446 if (trace & M_MATRIX) SAFE_ADD(M_count[i-1][j-1], count);
447 if (trace & Ix_MATRIX) SAFE_ADD(Ix_count[i-1][j-1], count);
448 if (trace & Iy_MATRIX) SAFE_ADD(Iy_count[i-1][j-1], count);
449 if (count == 0) count = 1; /* happens at M[0][0] only */
450 M_count[i][j] = count;
451 count = 0;
452 p = gaps[i][j].MIx;
453 if (p) {
454 while (1) {
455 gap = *p;
456 if (!gap) break;
457 SAFE_ADD(M_count[i-gap][j], count);
458 p++;
459 }
460 }
461 p = gaps[i][j].IyIx;
462 if (p) {
463 while (1) {
464 gap = *p;
465 if (!gap) break;
466 SAFE_ADD(Iy_count[i-gap][j], count);
467 p++;
468 }
469 }
470 Ix_count[i][j] = count;
471 count = 0;
472 p = gaps[i][j].MIy;
473 if (p) {
474 while (1) {
475 gap = *p;
476 if (!gap) break;
477 SAFE_ADD(M_count[i][j-gap], count);
478 p++;
479 }
480 }
481 p = gaps[i][j].IxIy;
482 if (p) {
483 while (1) {
484 gap = *p;
485 if (!gap) break;
486 SAFE_ADD(Ix_count[i][j-gap], count);
487 p++;
488 }
489 }
490 Iy_count[i][j] = count;
491 }
492 }
493 count = 0;
494 if (M[nA][nB].trace)
495 SAFE_ADD(M_count[nA][nB], count);
496 if (gaps[nA][nB].MIx[0] || gaps[nA][nB].IyIx[0])
497 SAFE_ADD(Ix_count[nA][nB], count);
498 if (gaps[nA][nB].MIy[0] || gaps[nA][nB].IxIy[0])
499 SAFE_ADD(Iy_count[nA][nB], count);
500 exit:
501 if (M_count) {
502 if (Ix_count) {
503 if (Iy_count) {
504 for (i = 0; i <= nA; i++) {
505 if (!M_count[i]) break;
506 PyMem_Free(M_count[i]);
507 if (!Ix_count[i]) break;
508 PyMem_Free(Ix_count[i]);
509 if (!Iy_count[i]) break;
510 PyMem_Free(Iy_count[i]);
511 }
512 PyMem_Free(Iy_count);
513 }
514 PyMem_Free(Ix_count);
515 }
516 PyMem_Free(M_count);
517 }
518 return count;
519 }
520
521 static Py_ssize_t
PathGenerator_waterman_smith_beyer_local_length(PathGenerator * self)522 PathGenerator_waterman_smith_beyer_local_length(PathGenerator* self)
523 {
524 int i;
525 int j;
526 int trace;
527 int* p;
528 int gap;
529 const int nA = self->nA;
530 const int nB = self->nB;
531 Trace** M = self->M;
532 TraceGapsWatermanSmithBeyer** gaps = self->gaps.waterman_smith_beyer;
533 Py_ssize_t term;
534 Py_ssize_t count = MEMORY_ERROR;
535 Py_ssize_t total = 0;
536 Py_ssize_t** M_count = NULL;
537 Py_ssize_t** Ix_count = NULL;
538 Py_ssize_t** Iy_count = NULL;
539 M_count = PyMem_Malloc((nA+1)*sizeof(Py_ssize_t*));
540 if (!M_count) goto exit;
541 Ix_count = PyMem_Malloc((nA+1)*sizeof(Py_ssize_t*));
542 if (!Ix_count) goto exit;
543 Iy_count = PyMem_Malloc((nA+1)*sizeof(Py_ssize_t*));
544 if (!Iy_count) goto exit;
545 for (i = 0; i <= nA; i++) {
546 M_count[i] = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
547 if (!M_count[i]) goto exit;
548 Ix_count[i] = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
549 if (!Ix_count[i]) goto exit;
550 Iy_count[i] = PyMem_Malloc((nB+1)*sizeof(Py_ssize_t));
551 if (!Iy_count[i]) goto exit;
552 }
553 for (i = 0; i <= nA; i++) {
554 for (j = 0; j <= nB; j++) {
555 count = 0;
556 trace = M[i][j].trace;
557 if (trace & M_MATRIX) SAFE_ADD(M_count[i-1][j-1], count);
558 if (trace & Ix_MATRIX) SAFE_ADD(Ix_count[i-1][j-1], count);
559 if (trace & Iy_MATRIX) SAFE_ADD(Iy_count[i-1][j-1], count);
560 if (count == 0 && (trace & STARTPOINT)) count = 1;
561 M_count[i][j] = count;
562 if (M[i][j].trace & ENDPOINT) SAFE_ADD(count, total);
563 count = 0;
564 p = gaps[i][j].MIx;
565 if (p) {
566 while (1) {
567 gap = *p;
568 if (!gap) break;
569 SAFE_ADD(M_count[i-gap][j], count);
570 p++;
571 }
572 }
573 p = gaps[i][j].IyIx;
574 if (p) {
575 while (1) {
576 gap = *p;
577 if (!gap) break;
578 SAFE_ADD(Iy_count[i-gap][j], count);
579 p++;
580 }
581 }
582 Ix_count[i][j] = count;
583 count = 0;
584 p = gaps[i][j].MIy;
585 if (p) {
586 while (1) {
587 gap = *p;
588 if (!gap) break;
589 SAFE_ADD(M_count[i][j-gap], count);
590 p++;
591 }
592 }
593 p = gaps[i][j].IxIy;
594 if (p) {
595 while (1) {
596 gap = *p;
597 if (!gap) break;
598 SAFE_ADD(Ix_count[i][j-gap], count);
599 p++;
600 }
601 }
602 Iy_count[i][j] = count;
603 }
604 }
605 count = total;
606 exit:
607 if (M_count) {
608 if (Ix_count) {
609 if (Iy_count) {
610 for (i = 0; i <= nA; i++) {
611 if (!M_count[i]) break;
612 PyMem_Free(M_count[i]);
613 if (!Ix_count[i]) break;
614 PyMem_Free(Ix_count[i]);
615 if (!Iy_count[i]) break;
616 PyMem_Free(Iy_count[i]);
617 }
618 PyMem_Free(Iy_count);
619 }
620 PyMem_Free(Ix_count);
621 }
622 PyMem_Free(M_count);
623 }
624 return count;
625 }
626
PathGenerator_length(PathGenerator * self)627 static Py_ssize_t PathGenerator_length(PathGenerator* self) {
628 Py_ssize_t length = self->length;
629 if (length == 0) {
630 switch (self->algorithm) {
631 case NeedlemanWunschSmithWaterman:
632 switch (self->mode) {
633 case Global:
634 length = PathGenerator_needlemanwunsch_length(self);
635 break;
636 case Local:
637 length = PathGenerator_smithwaterman_length(self);
638 break;
639 default:
640 /* should not happen, but some compilers complain that
641 * that length can be used uninitialized.
642 */
643 PyErr_SetString(PyExc_RuntimeError, "Unknown mode");
644 return -1;
645 }
646 break;
647 case Gotoh:
648 switch (self->mode) {
649 case Global:
650 length = PathGenerator_gotoh_global_length(self);
651 break;
652 case Local:
653 length = PathGenerator_gotoh_local_length(self);
654 break;
655 default:
656 /* should not happen, but some compilers complain that
657 * that length can be used uninitialized.
658 */
659 PyErr_SetString(PyExc_RuntimeError, "Unknown mode");
660 return -1;
661 }
662 break;
663 case WatermanSmithBeyer:
664 switch (self->mode) {
665 case Global:
666 length = PathGenerator_waterman_smith_beyer_global_length(self);
667 break;
668 case Local:
669 length = PathGenerator_waterman_smith_beyer_local_length(self);
670 break;
671 default:
672 /* should not happen, but some compilers complain that
673 * that length can be used uninitialized.
674 */
675 PyErr_SetString(PyExc_RuntimeError, "Unknown mode");
676 return -1;
677 }
678 break;
679 case Unknown:
680 default:
681 PyErr_SetString(PyExc_RuntimeError, "Unknown algorithm");
682 return -1;
683 }
684 self->length = length;
685 }
686 switch (length) {
687 case OVERFLOW_ERROR:
688 PyErr_Format(PyExc_OverflowError,
689 "number of optimal alignments is larger than %zd",
690 PY_SSIZE_T_MAX);
691 break;
692 case MEMORY_ERROR:
693 PyErr_SetNone(PyExc_MemoryError);
694 break;
695 default:
696 break;
697 }
698 return length;
699 }
700
701 static void
PathGenerator_dealloc(PathGenerator * self)702 PathGenerator_dealloc(PathGenerator* self)
703 {
704 int i;
705 const int nA = self->nA;
706 const Algorithm algorithm = self->algorithm;
707 Trace** M = self->M;
708 if (M) {
709 for (i = 0; i <= nA; i++) {
710 if (!M[i]) break;
711 PyMem_Free(M[i]);
712 }
713 PyMem_Free(M);
714 }
715 switch (algorithm) {
716 case NeedlemanWunschSmithWaterman:
717 break;
718 case Gotoh: {
719 TraceGapsGotoh** gaps = self->gaps.gotoh;
720 if (gaps) {
721 for (i = 0; i <= nA; i++) {
722 if (!gaps[i]) break;
723 PyMem_Free(gaps[i]);
724 }
725 PyMem_Free(gaps);
726 }
727 break;
728 }
729 case WatermanSmithBeyer: {
730 TraceGapsWatermanSmithBeyer** gaps = self->gaps.waterman_smith_beyer;
731 if (gaps) {
732 int j;
733 const int nB = self->nB;
734 int* trace;
735 for (i = 0; i <= nA; i++) {
736 if (!gaps[i]) break;
737 for (j = 0; j <= nB; j++) {
738 trace = gaps[i][j].MIx;
739 if (trace) PyMem_Free(trace);
740 trace = gaps[i][j].IyIx;
741 if (trace) PyMem_Free(trace);
742 trace = gaps[i][j].MIy;
743 if (trace) PyMem_Free(trace);
744 trace = gaps[i][j].IxIy;
745 if (trace) PyMem_Free(trace);
746 }
747 PyMem_Free(gaps[i]);
748 }
749 PyMem_Free(gaps);
750 }
751 break;
752 }
753 case Unknown:
754 default:
755 PyErr_WriteUnraisable((PyObject*)self);
756 break;
757 }
758 Py_TYPE(self)->tp_free((PyObject*)self);
759 }
760
PathGenerator_next_needlemanwunsch(PathGenerator * self)761 static PyObject* PathGenerator_next_needlemanwunsch(PathGenerator* self)
762 {
763 int i = 0;
764 int j = 0;
765 int path;
766 int trace = 0;
767 const int nA = self->nA;
768 const int nB = self->nB;
769 Trace** M = self->M;
770
771 path = M[i][j].path;
772 if (path == DONE) return NULL;
773 if (path == 0) {
774 /* Generate the first path. */
775 i = nA;
776 j = nB;
777 }
778 else {
779 /* We already have a path. Prune the path to see if there are
780 * any alternative paths. */
781 while (1) {
782 if (path == HORIZONTAL) {
783 trace = M[i][++j].trace;
784 if (trace & VERTICAL) {
785 M[--i][j].path = VERTICAL;
786 break;
787 }
788 if (trace & DIAGONAL) {
789 M[--i][--j].path = DIAGONAL;
790 break;
791 }
792 }
793 else if (path == VERTICAL) {
794 trace = M[++i][j].trace;
795 if (trace & DIAGONAL) {
796 M[--i][--j].path = DIAGONAL;
797 break;
798 }
799 }
800 else /* DIAGONAL */ {
801 i++;
802 j++;
803 }
804 path = M[i][j].path;
805 if (!path) {
806 /* we reached the end of the alignment without finding
807 * an alternative path */
808 M[0][0].path = DONE;
809 return NULL;
810 }
811 }
812 }
813 /* Follow the traceback until we reach the origin. */
814 while (1) {
815 trace = M[i][j].trace;
816 if (trace & HORIZONTAL) M[i][--j].path = HORIZONTAL;
817 else if (trace & VERTICAL) M[--i][j].path = VERTICAL;
818 else if (trace & DIAGONAL) M[--i][--j].path = DIAGONAL;
819 else break;
820 }
821 return PathGenerator_create_path(self, 0, 0);
822 }
823
PathGenerator_next_smithwaterman(PathGenerator * self)824 static PyObject* PathGenerator_next_smithwaterman(PathGenerator* self)
825 {
826 int trace = 0;
827 int i = self->iA;
828 int j = self->iB;
829 const int nA = self->nA;
830 const int nB = self->nB;
831 Trace** M = self->M;
832 int path = M[0][0].path;
833
834 if (path == DONE || path == NONE) return NULL;
835
836 path = M[i][j].path;
837 if (path) {
838 /* We already have a path. Prune the path to see if there are
839 * any alternative paths. */
840 while (1) {
841 if (path == HORIZONTAL) {
842 trace = M[i][++j].trace;
843 if (trace & VERTICAL) {
844 M[--i][j].path = VERTICAL;
845 break;
846 }
847 else if (trace & DIAGONAL) {
848 M[--i][--j].path = DIAGONAL;
849 break;
850 }
851 }
852 else if (path == VERTICAL) {
853 trace = M[++i][j].trace;
854 if (trace & DIAGONAL) {
855 M[--i][--j].path = DIAGONAL;
856 break;
857 }
858 }
859 else /* DIAGONAL */ {
860 i++;
861 j++;
862 }
863 path = M[i][j].path;
864 if (!path) break;
865 }
866 }
867
868 if (path) {
869 trace = M[i][j].trace;
870 } else {
871 /* Find a suitable end point for a path.
872 * Only allow end points ending at the M matrix. */
873 while (1) {
874 if (j < nB) j++;
875 else if (i < nA) {
876 i++;
877 j = 0;
878 }
879 else {
880 /* we reached the end of the sequences without finding
881 * an alternative path */
882 M[0][0].path = DONE;
883 return NULL;
884 }
885 trace = M[i][j].trace;
886 if (trace & ENDPOINT) {
887 trace &= DIAGONAL; /* exclude paths ending in a gap */
888 break;
889 }
890 }
891 M[i][j].path = 0;
892 }
893
894 /* Follow the traceback until we reach the origin. */
895 while (1) {
896 if (trace & HORIZONTAL) M[i][--j].path = HORIZONTAL;
897 else if (trace & VERTICAL) M[--i][j].path = VERTICAL;
898 else if (trace & DIAGONAL) M[--i][--j].path = DIAGONAL;
899 else if (trace & STARTPOINT) {
900 self->iA = i;
901 self->iB = j;
902 return PathGenerator_create_path(self, i, j);
903 }
904 else {
905 PyErr_SetString(PyExc_RuntimeError,
906 "Unexpected trace in PathGenerator_next_smithwaterman");
907 return NULL;
908 }
909 trace = M[i][j].trace;
910 }
911 }
912
PathGenerator_next_gotoh_global(PathGenerator * self)913 static PyObject* PathGenerator_next_gotoh_global(PathGenerator* self)
914 {
915 int i = 0;
916 int j = 0;
917 int m;
918 int path;
919 int trace = 0;
920 const int nA = self->nA;
921 const int nB = self->nB;
922 Trace** M = self->M;
923 TraceGapsGotoh** gaps = self->gaps.gotoh;
924
925 m = M_MATRIX;
926 path = M[i][j].path;
927 if (path == DONE) return NULL;
928 if (path == 0) {
929 i = nA;
930 j = nB;
931 }
932 else {
933 /* We already have a path. Prune the path to see if there are
934 * any alternative paths. */
935 while (1) {
936 path = M[i][j].path;
937 if (path == 0) {
938 switch (m) {
939 case M_MATRIX: m = Ix_MATRIX; break;
940 case Ix_MATRIX: m = Iy_MATRIX; break;
941 case Iy_MATRIX: m = 0; break;
942 }
943 break;
944 }
945 switch (path) {
946 case HORIZONTAL: trace = gaps[i][++j].Iy; break;
947 case VERTICAL: trace = gaps[++i][j].Ix; break;
948 case DIAGONAL: trace = M[++i][++j].trace; break;
949 }
950 switch (m) {
951 case M_MATRIX:
952 if (trace & Ix_MATRIX) {
953 m = Ix_MATRIX;
954 break;
955 }
956 case Ix_MATRIX:
957 if (trace & Iy_MATRIX) {
958 m = Iy_MATRIX;
959 break;
960 }
961 case Iy_MATRIX:
962 default:
963 switch (path) {
964 case HORIZONTAL: m = Iy_MATRIX; break;
965 case VERTICAL: m = Ix_MATRIX; break;
966 case DIAGONAL: m = M_MATRIX; break;
967 }
968 continue;
969 }
970 switch (path) {
971 case HORIZONTAL: j--; break;
972 case VERTICAL: i--; break;
973 case DIAGONAL: i--; j--; break;
974 }
975 M[i][j].path = path;
976 break;
977 }
978 }
979
980 if (path == 0) {
981 /* Generate a new path. */
982 switch (m) {
983 case M_MATRIX:
984 if (M[nA][nB].trace) {
985 /* m = M_MATRIX; */
986 break;
987 }
988 case Ix_MATRIX:
989 if (gaps[nA][nB].Ix) {
990 m = Ix_MATRIX;
991 break;
992 }
993 case Iy_MATRIX:
994 if (gaps[nA][nB].Iy) {
995 m = Iy_MATRIX;
996 break;
997 }
998 default:
999 /* exhausted this generator */
1000 M[0][0].path = DONE;
1001 return NULL;
1002 }
1003 }
1004
1005 switch (m) {
1006 case M_MATRIX:
1007 trace = M[i][j].trace;
1008 path = DIAGONAL;
1009 i--; j--;
1010 break;
1011 case Ix_MATRIX:
1012 trace = gaps[i][j].Ix;
1013 path = VERTICAL;
1014 i--;
1015 break;
1016 case Iy_MATRIX:
1017 trace = gaps[i][j].Iy;
1018 path = HORIZONTAL;
1019 j--;
1020 break;
1021 }
1022
1023 while (1) {
1024 if (trace & M_MATRIX) {
1025 trace = M[i][j].trace;
1026 M[i][j].path = path;
1027 path = DIAGONAL;
1028 i--; j--;
1029 }
1030 else if (trace & Ix_MATRIX) {
1031 M[i][j].path = path;
1032 trace = gaps[i][j].Ix;
1033 path = VERTICAL;
1034 i--;
1035 }
1036 else if (trace & Iy_MATRIX) {
1037 M[i][j].path = path;
1038 trace = gaps[i][j].Iy;
1039 path = HORIZONTAL;
1040 j--;
1041 }
1042 else break;
1043 }
1044 return PathGenerator_create_path(self, 0, 0);
1045 }
1046
PathGenerator_next_gotoh_local(PathGenerator * self)1047 static PyObject* PathGenerator_next_gotoh_local(PathGenerator* self)
1048 {
1049 int trace = 0;
1050 int i;
1051 int j;
1052 int m = M_MATRIX;
1053 int iA = self->iA;
1054 int iB = self->iB;
1055 const int nA = self->nA;
1056 const int nB = self->nB;
1057 Trace** M = self->M;
1058 TraceGapsGotoh** gaps = self->gaps.gotoh;
1059 int path = M[0][0].path;
1060
1061 if (path == DONE) return NULL;
1062
1063 path = M[iA][iB].path;
1064
1065 if (path) {
1066 i = iA;
1067 j = iB;
1068 while (1) {
1069 /* We already have a path. Prune the path to see if there are
1070 * any alternative paths. */
1071 path = M[i][j].path;
1072 if (path == 0) {
1073 m = M_MATRIX;
1074 iA = i;
1075 iB = j;
1076 break;
1077 }
1078 switch (path) {
1079 case HORIZONTAL: trace = gaps[i][++j].Iy; break;
1080 case VERTICAL: trace = gaps[++i][j].Ix; break;
1081 case DIAGONAL: trace = M[++i][++j].trace; break;
1082 }
1083 switch (m) {
1084 case M_MATRIX:
1085 if (trace & Ix_MATRIX) {
1086 m = Ix_MATRIX;
1087 break;
1088 }
1089 case Ix_MATRIX:
1090 if (trace & Iy_MATRIX) {
1091 m = Iy_MATRIX;
1092 break;
1093 }
1094 case Iy_MATRIX:
1095 default:
1096 switch (path) {
1097 case HORIZONTAL: m = Iy_MATRIX; break;
1098 case VERTICAL: m = Ix_MATRIX; break;
1099 case DIAGONAL: m = M_MATRIX; break;
1100 }
1101 continue;
1102 }
1103 switch (path) {
1104 case HORIZONTAL: j--; break;
1105 case VERTICAL: i--; break;
1106 case DIAGONAL: i--; j--; break;
1107 }
1108 M[i][j].path = path;
1109 break;
1110 }
1111 }
1112
1113 if (path == 0) {
1114 /* Find the end point for a new path. */
1115 while (1) {
1116 if (iB < nB) iB++;
1117 else if (iA < nA) {
1118 iA++;
1119 iB = 0;
1120 }
1121 else {
1122 /* we reached the end of the alignment without finding
1123 * an alternative path */
1124 M[0][0].path = DONE;
1125 return NULL;
1126 }
1127 if (M[iA][iB].trace & ENDPOINT) {
1128 M[iA][iB].path = 0;
1129 break;
1130 }
1131 }
1132 m = M_MATRIX;
1133 i = iA;
1134 j = iB;
1135 }
1136
1137 while (1) {
1138 switch (m) {
1139 case M_MATRIX: trace = M[i][j].trace; break;
1140 case Ix_MATRIX: trace = gaps[i][j].Ix; break;
1141 case Iy_MATRIX: trace = gaps[i][j].Iy; break;
1142 }
1143 if (trace == STARTPOINT) {
1144 self->iA = i;
1145 self->iB = j;
1146 return PathGenerator_create_path(self, i, j);
1147 }
1148 switch (m) {
1149 case M_MATRIX:
1150 path = DIAGONAL;
1151 i--;
1152 j--;
1153 break;
1154 case Ix_MATRIX:
1155 path = VERTICAL;
1156 i--;
1157 break;
1158 case Iy_MATRIX:
1159 path = HORIZONTAL;
1160 j--;
1161 break;
1162 }
1163 if (trace & M_MATRIX) m = M_MATRIX;
1164 else if (trace & Ix_MATRIX) m = Ix_MATRIX;
1165 else if (trace & Iy_MATRIX) m = Iy_MATRIX;
1166 else {
1167 PyErr_SetString(PyExc_RuntimeError,
1168 "Unexpected trace in PathGenerator_next_gotoh_local");
1169 return NULL;
1170 }
1171 M[i][j].path = path;
1172 }
1173 return NULL;
1174 }
1175
1176 static PyObject*
PathGenerator_next_waterman_smith_beyer_global(PathGenerator * self)1177 PathGenerator_next_waterman_smith_beyer_global(PathGenerator* self)
1178 {
1179 int i = 0, j = 0;
1180 int iA, iB;
1181 int trace;
1182 int* gapM;
1183 int* gapXY;
1184
1185 int m = M_MATRIX;
1186 const int nA = self->nA;
1187 const int nB = self->nB;
1188 Trace** M = self->M;
1189 TraceGapsWatermanSmithBeyer** gaps = self->gaps.waterman_smith_beyer;
1190
1191 int gap;
1192 int path = M[0][0].path;
1193
1194 if (path == DONE) return NULL;
1195
1196 if (path) {
1197 /* We already have a path. Prune the path to see if there are
1198 * any alternative paths. */
1199 while (1) {
1200 if (!path) {
1201 m <<= 1;
1202 break;
1203 }
1204 switch (path) {
1205 case HORIZONTAL:
1206 iA = i;
1207 iB = j;
1208 while (M[i][iB].path == HORIZONTAL) iB++;
1209 break;
1210 case VERTICAL:
1211 iA = i;
1212 while (M[iA][j].path == VERTICAL) iA++;
1213 iB = j;
1214 break;
1215 case DIAGONAL:
1216 iA = i + 1;
1217 iB = j + 1;
1218 break;
1219 default:
1220 PyErr_SetString(PyExc_RuntimeError,
1221 "Unexpected path in PathGenerator_next_waterman_smith_beyer_global");
1222 return NULL;
1223 }
1224 if (i == iA) { /* HORIZONTAL */
1225 gapM = gaps[iA][iB].MIy;
1226 gapXY = gaps[iA][iB].IxIy;
1227 if (m == M_MATRIX) {
1228 gap = iB - j;
1229 while (*gapM != gap) gapM++;
1230 gapM++;
1231 gap = *gapM;
1232 if (gap) {
1233 j = iB - gap;
1234 while (j < iB) M[i][--iB].path = HORIZONTAL;
1235 break;
1236 }
1237 } else if (m == Ix_MATRIX) {
1238 gap = iB - j;
1239 while (*gapXY != gap) gapXY++;
1240 gapXY++;
1241 }
1242 gap = *gapXY;
1243 if (gap) {
1244 m = Ix_MATRIX;
1245 j = iB - gap;
1246 while (j < iB) M[i][--iB].path = HORIZONTAL;
1247 break;
1248 }
1249 /* no alternative found; continue pruning */
1250 m = Iy_MATRIX;
1251 j = iB;
1252 }
1253 else if (j == iB) { /* VERTICAL */
1254 gapM = gaps[iA][iB].MIx;
1255 gapXY = gaps[iA][iB].IyIx;
1256 if (m == M_MATRIX) {
1257 gap = iA - i;
1258 while (*gapM != gap) gapM++;
1259 gapM++;
1260 gap = *gapM;
1261 if (gap) {
1262 i = iA - gap;
1263 while (i < iA) M[--iA][j].path = VERTICAL;
1264 break;
1265 }
1266 } else if (m == Iy_MATRIX) {
1267 gap = iA - i;
1268 while (*gapXY != gap) gapXY++;
1269 gapXY++;
1270 }
1271 gap = *gapXY;
1272 if (gap) {
1273 m = Iy_MATRIX;
1274 i = iA - gap;
1275 while (i < iA) M[--iA][j].path = VERTICAL;
1276 break;
1277 }
1278 /* no alternative found; continue pruning */
1279 m = Ix_MATRIX;
1280 i = iA;
1281 }
1282 else { /* DIAGONAL */
1283 i = iA - 1;
1284 j = iB - 1;
1285 trace = M[iA][iB].trace;
1286 switch (m) {
1287 case M_MATRIX:
1288 if (trace & Ix_MATRIX) {
1289 m = Ix_MATRIX;
1290 M[i][j].path = DIAGONAL;
1291 break;
1292 }
1293 case Ix_MATRIX:
1294 if (trace & Iy_MATRIX) {
1295 m = Iy_MATRIX;
1296 M[i][j].path = DIAGONAL;
1297 break;
1298 }
1299 case Iy_MATRIX:
1300 default:
1301 /* no alternative found; continue pruning */
1302 m = M_MATRIX;
1303 i = iA;
1304 j = iB;
1305 path = M[i][j].path;
1306 continue;
1307 }
1308 /* alternative found; build path until starting point */
1309 break;
1310 }
1311 path = M[i][j].path;
1312 }
1313 }
1314
1315 if (!path) {
1316 /* Find a suitable end point for a path. */
1317 switch (m) {
1318 case M_MATRIX:
1319 if (M[nA][nB].trace) {
1320 /* m = M_MATRIX; */
1321 break;
1322 }
1323 case Ix_MATRIX:
1324 if (gaps[nA][nB].MIx[0] || gaps[nA][nB].IyIx[0]) {
1325 m = Ix_MATRIX;
1326 break;
1327 }
1328 case Iy_MATRIX:
1329 if (gaps[nA][nB].MIy[0] || gaps[nA][nB].IxIy[0]) {
1330 m = Iy_MATRIX;
1331 break;
1332 }
1333 default:
1334 M[0][0].path = DONE;
1335 return NULL;
1336 }
1337 i = nA;
1338 j = nB;
1339 }
1340
1341 /* Follow the traceback until we reach the origin. */
1342 while (1) {
1343 switch (m) {
1344 case M_MATRIX:
1345 trace = M[i][j].trace;
1346 if (trace & M_MATRIX) m = M_MATRIX;
1347 else if (trace & Ix_MATRIX) m = Ix_MATRIX;
1348 else if (trace & Iy_MATRIX) m = Iy_MATRIX;
1349 else return PathGenerator_create_path(self, i, j);
1350 i--;
1351 j--;
1352 M[i][j].path = DIAGONAL;
1353 break;
1354 case Ix_MATRIX:
1355 gap = gaps[i][j].MIx[0];
1356 if (gap) m = M_MATRIX;
1357 else {
1358 gap = gaps[i][j].IyIx[0];
1359 m = Iy_MATRIX;
1360 }
1361 iA = i - gap;
1362 while (iA < i) M[--i][j].path = VERTICAL;
1363 M[i][j].path = VERTICAL;
1364 break;
1365 case Iy_MATRIX:
1366 gap = gaps[i][j].MIy[0];
1367 if (gap) m = M_MATRIX;
1368 else {
1369 gap = gaps[i][j].IxIy[0];
1370 m = Ix_MATRIX;
1371 }
1372 iB = j - gap;
1373 while (iB < j) M[i][--j].path = HORIZONTAL;
1374 M[i][j].path = HORIZONTAL;
1375 break;
1376 }
1377 }
1378 }
1379
1380 static PyObject*
PathGenerator_next_waterman_smith_beyer_local(PathGenerator * self)1381 PathGenerator_next_waterman_smith_beyer_local(PathGenerator* self)
1382 {
1383 int i, j, m;
1384 int trace = 0;
1385 int* gapM;
1386 int* gapXY;
1387
1388 int iA = self->iA;
1389 int iB = self->iB;
1390 const int nA = self->nA;
1391 const int nB = self->nB;
1392 Trace** M = self->M;
1393 TraceGapsWatermanSmithBeyer** gaps = self->gaps.waterman_smith_beyer;
1394
1395 int gap;
1396 int path = M[0][0].path;
1397
1398 if (path == DONE) return NULL;
1399 m = 0;
1400 path = M[iA][iB].path;
1401 if (path) {
1402 /* We already have a path. Prune the path to see if there are
1403 * any alternative paths. */
1404 m = M_MATRIX;
1405 i = iA;
1406 j = iB;
1407 while (1) {
1408 path = M[i][j].path;
1409 switch (path) {
1410 case HORIZONTAL:
1411 iA = i;
1412 iB = j;
1413 while (M[i][iB].path == HORIZONTAL) iB++;
1414 break;
1415 case VERTICAL:
1416 iA = i;
1417 iB = j;
1418 while (M[iA][j].path == VERTICAL) iA++;
1419 break;
1420 case DIAGONAL:
1421 iA = i + 1;
1422 iB = j + 1;
1423 break;
1424 default:
1425 iA = -1;
1426 break;
1427 }
1428 if (iA < 0) {
1429 m = 0;
1430 iA = i;
1431 iB = j;
1432 break;
1433 }
1434 if (i == iA) { /* HORIZONTAL */
1435 gapM = gaps[iA][iB].MIy;
1436 gapXY = gaps[iA][iB].IxIy;
1437 if (m == M_MATRIX) {
1438 gap = iB - j;
1439 while (*gapM != gap) gapM++;
1440 gapM++;
1441 gap = *gapM;
1442 if (gap) {
1443 j = iB - gap;
1444 while (j < iB) M[i][--iB].path = HORIZONTAL;
1445 break;
1446 }
1447 } else if (m == Ix_MATRIX) {
1448 gap = iB - j;
1449 while (*gapXY != gap) gapXY++;
1450 gapXY++;
1451 }
1452 gap = *gapXY;
1453 if (gap) {
1454 m = Ix_MATRIX;
1455 j = iB - gap;
1456 M[i][j].path = HORIZONTAL;
1457 while (iB > j) M[i][--iB].path = HORIZONTAL;
1458 break;
1459 }
1460 /* no alternative found; continue pruning */
1461 m = Iy_MATRIX;
1462 j = iB;
1463 }
1464 else if (j == iB) { /* VERTICAL */
1465 gapM = gaps[iA][iB].MIx;
1466 gapXY = gaps[iA][iB].IyIx;
1467 if (m == M_MATRIX) {
1468 gap = iA - i;
1469 while (*gapM != gap) gapM++;
1470 gapM++;
1471 gap = *gapM;
1472 if (gap) {
1473 i = iA - gap;
1474 while (i < iA) M[--iA][j].path = VERTICAL;
1475 break;
1476 }
1477 } else if (m == Iy_MATRIX) {
1478 gap = iA - i;
1479 while (*gapXY != gap) gapXY++;
1480 gapXY++;
1481 }
1482 gap = *gapXY;
1483 if (gap) {
1484 m = Iy_MATRIX;
1485 i = iA - gap;
1486 M[i][j].path = VERTICAL;
1487 while (iA > i) M[--iA][j].path = VERTICAL;
1488 break;
1489 }
1490 /* no alternative found; continue pruning */
1491 m = Ix_MATRIX;
1492 i = iA;
1493 }
1494 else { /* DIAGONAL */
1495 i = iA - 1;
1496 j = iB - 1;
1497 trace = M[iA][iB].trace;
1498 switch (m) {
1499 case M_MATRIX:
1500 if (trace & Ix_MATRIX) {
1501 m = Ix_MATRIX;
1502 M[i][j].path = DIAGONAL;
1503 break;
1504 }
1505 case Ix_MATRIX:
1506 if (trace & Iy_MATRIX) {
1507 m = Iy_MATRIX;
1508 M[i][j].path = DIAGONAL;
1509 break;
1510 }
1511 case Iy_MATRIX:
1512 default:
1513 /* no alternative found; continue pruning */
1514 m = M_MATRIX;
1515 i = iA;
1516 j = iB;
1517 continue;
1518 }
1519 /* alternative found; build path until starting point */
1520 break;
1521 }
1522 }
1523 }
1524
1525 if (m == 0) {
1526 /* We are at [nA][nB]. Find a suitable end point for a path. */
1527 while (1) {
1528 if (iB < nB) iB++;
1529 else if (iA < nA) {
1530 iA++;
1531 iB = 0;
1532 }
1533 else {
1534 /* exhausted this generator */
1535 M[0][0].path = DONE;
1536 return NULL;
1537 }
1538 if (M[iA][iB].trace & ENDPOINT) break;
1539 }
1540 M[iA][iB].path = 0;
1541 m = M_MATRIX;
1542 i = iA;
1543 j = iB;
1544 }
1545
1546 /* Follow the traceback until we reach the origin. */
1547 while (1) {
1548 switch (m) {
1549 case Ix_MATRIX:
1550 gapM = gaps[i][j].MIx;
1551 gapXY = gaps[i][j].IyIx;
1552 iB = j;
1553 gap = *gapM;
1554 if (gap) m = M_MATRIX;
1555 else {
1556 gap = *gapXY;
1557 m = Iy_MATRIX;
1558 }
1559 iA = i - gap;
1560 while (i > iA) M[--i][iB].path = VERTICAL;
1561 break;
1562 case Iy_MATRIX:
1563 gapM = gaps[i][j].MIy;
1564 gapXY = gaps[i][j].IxIy;
1565 iA = i;
1566 gap = *gapM;
1567 if (gap) m = M_MATRIX;
1568 else {
1569 gap = *gapXY;
1570 m = Ix_MATRIX;
1571 }
1572 iB = j - gap;
1573 while (j > iB) M[iA][--j].path = HORIZONTAL;
1574 break;
1575 case M_MATRIX:
1576 iA = i-1;
1577 iB = j-1;
1578 trace = M[i][j].trace;
1579 if (trace & M_MATRIX) m = M_MATRIX;
1580 else if (trace & Ix_MATRIX) m = Ix_MATRIX;
1581 else if (trace & Iy_MATRIX) m = Iy_MATRIX;
1582 else if (trace == STARTPOINT) {
1583 self->iA = i;
1584 self->iB = j;
1585 return PathGenerator_create_path(self, i, j);
1586 }
1587 else {
1588 PyErr_SetString(PyExc_RuntimeError,
1589 "Unexpected trace in PathGenerator_next_waterman_smith_beyer_local");
1590 return NULL;
1591 }
1592 M[iA][iB].path = DIAGONAL;
1593 break;
1594 }
1595 i = iA;
1596 j = iB;
1597 }
1598 }
1599
1600 static PyObject *
PathGenerator_next(PathGenerator * self)1601 PathGenerator_next(PathGenerator* self)
1602 {
1603 const Mode mode = self->mode;
1604 const Algorithm algorithm = self->algorithm;
1605 switch (algorithm) {
1606 case NeedlemanWunschSmithWaterman:
1607 switch (mode) {
1608 case Global:
1609 return PathGenerator_next_needlemanwunsch(self);
1610 case Local:
1611 return PathGenerator_next_smithwaterman(self);
1612 }
1613 case Gotoh:
1614 switch (mode) {
1615 case Global:
1616 return PathGenerator_next_gotoh_global(self);
1617 case Local:
1618 return PathGenerator_next_gotoh_local(self);
1619 }
1620 case WatermanSmithBeyer:
1621 switch (mode) {
1622 case Global:
1623 return PathGenerator_next_waterman_smith_beyer_global(self);
1624 case Local:
1625 return PathGenerator_next_waterman_smith_beyer_local(self);
1626 }
1627 case Unknown:
1628 default:
1629 PyErr_SetString(PyExc_RuntimeError, "Unknown algorithm");
1630 return NULL;
1631 }
1632 }
1633
1634 static const char PathGenerator_reset__doc__[] = "reset the iterator";
1635
1636 static PyObject*
PathGenerator_reset(PathGenerator * self)1637 PathGenerator_reset(PathGenerator* self)
1638 {
1639 switch (self->mode) {
1640 case Local:
1641 self->iA = 0;
1642 self->iB = 0;
1643 case Global: {
1644 Trace** M = self->M;
1645 switch (self->algorithm) {
1646 case NeedlemanWunschSmithWaterman:
1647 case Gotoh: {
1648 if (M[0][0].path != NONE) M[0][0].path = 0;
1649 break;
1650 }
1651 case WatermanSmithBeyer: {
1652 M[0][0].path = 0;
1653 break;
1654 }
1655 case Unknown:
1656 default:
1657 break;
1658 }
1659 }
1660 }
1661 Py_INCREF(Py_None);
1662 return Py_None;
1663 }
1664
1665 static PyMethodDef PathGenerator_methods[] = {
1666 {"reset",
1667 (PyCFunction)PathGenerator_reset,
1668 METH_NOARGS,
1669 PathGenerator_reset__doc__
1670 },
1671 {NULL} /* Sentinel */
1672 };
1673
1674 static PySequenceMethods PathGenerator_as_sequence = {
1675 (lenfunc)PathGenerator_length, /* sq_length */
1676 NULL, /* sq_concat */
1677 NULL, /* sq_repeat */
1678 NULL, /* sq_item */
1679 NULL, /* sq_ass_item */
1680 NULL, /* sq_contains */
1681 NULL, /* sq_inplace_concat */
1682 NULL, /* sq_inplace_repeat */
1683 };
1684
1685 static PyTypeObject PathGenerator_Type = {
1686 PyVarObject_HEAD_INIT(NULL, 0)
1687 "Path generator", /* tp_name */
1688 sizeof(PathGenerator), /* tp_basicsize */
1689 0, /* tp_itemsize */
1690 (destructor)PathGenerator_dealloc, /* tp_dealloc */
1691 0, /* tp_print */
1692 0, /* tp_getattr */
1693 0, /* tp_setattr */
1694 0, /* tp_reserved */
1695 0, /* tp_repr */
1696 0, /* tp_as_number */
1697 &PathGenerator_as_sequence, /* tp_as_sequence */
1698 0, /* tp_as_mapping */
1699 0, /* tp_hash */
1700 0, /* tp_call */
1701 0, /* tp_str */
1702 0, /* tp_getattro */
1703 0, /* tp_setattro */
1704 0, /* tp_as_buffer */
1705 Py_TPFLAGS_DEFAULT, /* tp_flags */
1706 0, /* tp_doc */
1707 0, /* tp_traverse */
1708 0, /* tp_clear */
1709 0, /* tp_richcompare */
1710 0, /* tp_weaklistoffset */
1711 PyObject_SelfIter, /* tp_iter */
1712 (iternextfunc)PathGenerator_next, /* tp_iternext */
1713 PathGenerator_methods, /* tp_methods */
1714 };
1715
1716 typedef struct {
1717 PyObject_HEAD
1718 Mode mode;
1719 Algorithm algorithm;
1720 double match;
1721 double mismatch;
1722 double epsilon;
1723 double target_internal_open_gap_score;
1724 double target_internal_extend_gap_score;
1725 double target_left_open_gap_score;
1726 double target_left_extend_gap_score;
1727 double target_right_open_gap_score;
1728 double target_right_extend_gap_score;
1729 double query_internal_open_gap_score;
1730 double query_internal_extend_gap_score;
1731 double query_left_open_gap_score;
1732 double query_left_extend_gap_score;
1733 double query_right_open_gap_score;
1734 double query_right_extend_gap_score;
1735 PyObject* target_gap_function;
1736 PyObject* query_gap_function;
1737 Py_buffer substitution_matrix;
1738 PyObject* alphabet;
1739 int* mapping;
1740 int wildcard;
1741 } Aligner;
1742
1743
1744 static Py_ssize_t
set_alphabet(Aligner * self,PyObject * alphabet)1745 set_alphabet(Aligner* self, PyObject* alphabet)
1746 {
1747 Py_ssize_t size;
1748 if (alphabet == Py_None) {
1749 if (self->alphabet) {
1750 Py_DECREF(self->alphabet);
1751 self->alphabet = NULL;
1752 }
1753 if (self->mapping) {
1754 PyMem_Free(self->mapping);
1755 self->mapping = NULL;
1756 }
1757 return 0;
1758 }
1759 else if (PyUnicode_Check(alphabet)) {
1760 int* mapping;
1761 int i;
1762 int n;
1763 int kind;
1764 void* characters;
1765 if (PyUnicode_READY(alphabet) == -1) return -1;
1766 size = PyUnicode_GET_LENGTH(alphabet);
1767 if (size == 0) {
1768 PyErr_SetString(PyExc_ValueError, "alphabet has zero length");
1769 return -1;
1770 }
1771 kind = PyUnicode_KIND(alphabet);
1772 switch (kind) {
1773 case PyUnicode_1BYTE_KIND: {
1774 n = 1 << 8 * sizeof(Py_UCS1);
1775 break;
1776 }
1777 case PyUnicode_2BYTE_KIND: {
1778 n = 1 << 8 * sizeof(Py_UCS2);
1779 break;
1780 }
1781 case PyUnicode_4BYTE_KIND: {
1782 n = 0x110000; /* Maximum code point in Unicode 6.0
1783 * is 0x10ffff = 1114111 */
1784 break;
1785 }
1786 case PyUnicode_WCHAR_KIND:
1787 default:
1788 PyErr_SetString(PyExc_ValueError, "could not interpret alphabet");
1789 return -1;
1790 }
1791 characters = PyUnicode_DATA(alphabet);
1792 mapping = PyMem_Malloc(n*sizeof(int));
1793 if (!mapping) return -1;
1794 for (i = 0; i < n; i++) mapping[i] = MISSING_LETTER;
1795 for (i = 0; i < size; i++) {
1796 Py_UCS4 character = PyUnicode_READ(kind, characters, i);
1797 if (mapping[character] != MISSING_LETTER) {
1798 PyObject* c = PyUnicode_FromKindAndData(kind, &character, 1);
1799 PyErr_Format(PyExc_ValueError,
1800 "alphabet contains '%S' more than once", c);
1801 Py_XDECREF(c);
1802 PyMem_Free(mapping);
1803 return -1;
1804 }
1805 mapping[character] = i;
1806 }
1807 Py_INCREF(alphabet);
1808 if (self->mapping) PyMem_Free(self->mapping);
1809 self->mapping = mapping;
1810 }
1811 else {
1812 /* alphabet is not a string; cannot use mapping */
1813 PyObject* sequence = PySequence_Fast(alphabet,
1814 "alphabet should support the sequence protocol (e.g.,\n"
1815 "strings, lists, and tuples can be valid alphabets).");
1816 if (!sequence) return -1;
1817 size = PySequence_Fast_GET_SIZE(sequence);
1818 Py_DECREF(sequence);
1819 if (self->mapping) {
1820 PyMem_Free(self->mapping);
1821 self->mapping = NULL;
1822 }
1823 Py_INCREF(alphabet);
1824 }
1825 Py_XDECREF(self->alphabet);
1826 self->alphabet = alphabet;
1827 return size;
1828 }
1829
1830 static int
Aligner_init(Aligner * self,PyObject * args,PyObject * kwds)1831 Aligner_init(Aligner *self, PyObject *args, PyObject *kwds)
1832 {
1833 self->mode = Global;
1834 self->match = 1.0;
1835 self->mismatch = 0.0;
1836 self->epsilon = 1.e-6;
1837 self->target_internal_open_gap_score = 0;
1838 self->target_internal_extend_gap_score = 0;
1839 self->query_internal_open_gap_score = 0;
1840 self->query_internal_extend_gap_score = 0;
1841 self->target_left_open_gap_score = 0;
1842 self->target_left_extend_gap_score = 0;
1843 self->target_right_open_gap_score = 0;
1844 self->target_right_extend_gap_score = 0;
1845 self->query_left_open_gap_score = 0;
1846 self->query_left_extend_gap_score = 0;
1847 self->query_right_open_gap_score = 0;
1848 self->query_right_extend_gap_score = 0;
1849 self->target_gap_function = NULL;
1850 self->query_gap_function = NULL;
1851 self->substitution_matrix.obj = NULL;
1852 self->substitution_matrix.buf = NULL;
1853 self->algorithm = Unknown;
1854 self->alphabet = NULL;
1855 self->mapping = NULL;
1856 self->wildcard = -1;
1857 return 0;
1858 }
1859
1860 static void
Aligner_dealloc(Aligner * self)1861 Aligner_dealloc(Aligner* self)
1862 { Py_XDECREF(self->target_gap_function);
1863 Py_XDECREF(self->query_gap_function);
1864 if (self->substitution_matrix.obj) PyBuffer_Release(&self->substitution_matrix);
1865 Py_XDECREF(self->alphabet);
1866 Py_XDECREF(self->mapping);
1867 Py_TYPE(self)->tp_free((PyObject*)self);
1868 }
1869
1870 static PyObject*
Aligner_repr(Aligner * self)1871 Aligner_repr(Aligner* self)
1872 {
1873 const char text[] = "Pairwise aligner, implementing the Needleman-Wunsch, Smith-Waterman, Gotoh, and Waterman-Smith-Beyer global and local alignment algorithms";
1874 return PyUnicode_FromString(text);
1875 }
1876
1877 static PyObject*
Aligner_str(Aligner * self)1878 Aligner_str(Aligner* self)
1879 {
1880 char text[1024];
1881 char* p = text;
1882 PyObject* substitution_matrix = self->substitution_matrix.obj;
1883 void* args[3];
1884 int n = 0;
1885 PyObject* wildcard = NULL;
1886 PyObject* s;
1887
1888 p += sprintf(p, "Pairwise sequence aligner with parameters\n");
1889 if (substitution_matrix) {
1890 p += sprintf(p, " substitution_matrix: <%s object at %p>\n",
1891 Py_TYPE(substitution_matrix)->tp_name,
1892 substitution_matrix);
1893 } else {
1894 if (self->wildcard == -1) {
1895 p += sprintf(p, " wildcard: None\n");
1896 }
1897 else {
1898 wildcard = PyUnicode_FromKindAndData(PyUnicode_4BYTE_KIND,
1899 &self->wildcard, 1);
1900 if (!wildcard) return NULL;
1901 p += sprintf(p, " wildcard: '%%U'\n");
1902 args[n++] = wildcard;
1903 }
1904 p += sprintf(p, " match_score: %f\n", self->match);
1905 p += sprintf(p, " mismatch_score: %f\n", self->mismatch);
1906 }
1907 if (self->target_gap_function) {
1908 p += sprintf(p, " target_gap_function: %%R\n");
1909 args[n++] = self->target_gap_function;
1910 }
1911 else {
1912 p += sprintf(p, " target_internal_open_gap_score: %f\n",
1913 self->target_internal_open_gap_score);
1914 p += sprintf(p, " target_internal_extend_gap_score: %f\n",
1915 self->target_internal_extend_gap_score);
1916 p += sprintf(p, " target_left_open_gap_score: %f\n",
1917 self->target_left_open_gap_score);
1918 p += sprintf(p, " target_left_extend_gap_score: %f\n",
1919 self->target_left_extend_gap_score);
1920 p += sprintf(p, " target_right_open_gap_score: %f\n",
1921 self->target_right_open_gap_score);
1922 p += sprintf(p, " target_right_extend_gap_score: %f\n",
1923 self->target_right_extend_gap_score);
1924 }
1925 if (self->query_gap_function) {
1926 p += sprintf(p, " query_gap_function: %%R\n");
1927 args[n++] = self->query_gap_function;
1928 }
1929 else {
1930 p += sprintf(p, " query_internal_open_gap_score: %f\n",
1931 self->query_internal_open_gap_score);
1932 p += sprintf(p, " query_internal_extend_gap_score: %f\n",
1933 self->query_internal_extend_gap_score);
1934 p += sprintf(p, " query_left_open_gap_score: %f\n",
1935 self->query_left_open_gap_score);
1936 p += sprintf(p, " query_left_extend_gap_score: %f\n",
1937 self->query_left_extend_gap_score);
1938 p += sprintf(p, " query_right_open_gap_score: %f\n",
1939 self->query_right_open_gap_score);
1940 p += sprintf(p, " query_right_extend_gap_score: %f\n",
1941 self->query_right_extend_gap_score);
1942 }
1943 switch (self->mode) {
1944 case Global: sprintf(p, " mode: global\n"); break;
1945 case Local: sprintf(p, " mode: local\n"); break;
1946 }
1947 s = PyUnicode_FromFormat(text, args[0], args[1], args[2]);
1948 Py_XDECREF(wildcard);
1949 return s;
1950 }
1951
1952 static char Aligner_mode__doc__[] = "alignment mode ('global' or 'local')";
1953
1954 static PyObject*
Aligner_get_mode(Aligner * self,void * closure)1955 Aligner_get_mode(Aligner* self, void* closure)
1956 { const char* message = NULL;
1957 switch (self->mode) {
1958 case Global: message = "global"; break;
1959 case Local: message = "local"; break;
1960 }
1961 return PyUnicode_FromString(message);
1962 }
1963
1964 static int
Aligner_set_mode(Aligner * self,PyObject * value,void * closure)1965 Aligner_set_mode(Aligner* self, PyObject* value, void* closure)
1966 {
1967 if (PyUnicode_Check(value)) {
1968 if (PyUnicode_CompareWithASCIIString(value, "global") == 0) {
1969 self->mode = Global;
1970 return 0;
1971 }
1972 if (PyUnicode_CompareWithASCIIString(value, "local") == 0) {
1973 self->mode = Local;
1974 return 0;
1975 }
1976 }
1977 PyErr_SetString(PyExc_ValueError,
1978 "invalid mode (expected 'global' or 'local'");
1979 return -1;
1980 }
1981
1982 static char Aligner_match_score__doc__[] = "match score";
1983
1984 static PyObject*
Aligner_get_match_score(Aligner * self,void * closure)1985 Aligner_get_match_score(Aligner* self, void* closure)
1986 { if (self->substitution_matrix.obj) {
1987 Py_INCREF(Py_None);
1988 return Py_None;
1989 }
1990 return PyFloat_FromDouble(self->match);
1991 }
1992
1993 static int
Aligner_set_match_score(Aligner * self,PyObject * value,void * closure)1994 Aligner_set_match_score(Aligner* self, PyObject* value, void* closure)
1995 {
1996 const double match = PyFloat_AsDouble(value);
1997 if (PyErr_Occurred()) {
1998 PyErr_SetString(PyExc_ValueError, "invalid match score");
1999 return -1;
2000 }
2001 if (self->substitution_matrix.obj) {
2002 if (set_alphabet(self, Py_None) < 0) return -1;
2003 PyBuffer_Release(&self->substitution_matrix);
2004 }
2005 self->match = match;
2006 return 0;
2007 }
2008
2009 static char Aligner_mismatch_score__doc__[] = "mismatch score";
2010
2011 static PyObject*
Aligner_get_mismatch_score(Aligner * self,void * closure)2012 Aligner_get_mismatch_score(Aligner* self, void* closure)
2013 { if (self->substitution_matrix.obj) {
2014 Py_INCREF(Py_None);
2015 return Py_None;
2016 }
2017 return PyFloat_FromDouble(self->mismatch);
2018 }
2019
2020 static int
Aligner_set_mismatch_score(Aligner * self,PyObject * value,void * closure)2021 Aligner_set_mismatch_score(Aligner* self, PyObject* value, void* closure)
2022 {
2023 const double mismatch = PyFloat_AsDouble(value);
2024 if (PyErr_Occurred()) {
2025 PyErr_SetString(PyExc_ValueError, "invalid mismatch score");
2026 return -1;
2027 }
2028 if (self->substitution_matrix.obj) {
2029 if (set_alphabet(self, Py_None) < 0) return -1;
2030 PyBuffer_Release(&self->substitution_matrix);
2031 }
2032 self->mismatch = mismatch;
2033 return 0;
2034 }
2035
2036 static char Aligner_substitution_matrix__doc__[] = "substitution_matrix";
2037
2038 static PyObject*
Aligner_get_substitution_matrix(Aligner * self,void * closure)2039 Aligner_get_substitution_matrix(Aligner* self, void* closure)
2040 { PyObject* object = self->substitution_matrix.obj;
2041 if (!object) object = Py_None;
2042 Py_INCREF(object);
2043 return object;
2044 }
2045
2046 static int
Aligner_set_substitution_matrix(Aligner * self,PyObject * values,void * closure)2047 Aligner_set_substitution_matrix(Aligner* self, PyObject* values, void* closure)
2048 {
2049 PyObject* alphabet;
2050 Py_ssize_t size = -1;
2051 Py_buffer view;
2052 const int flag = PyBUF_FORMAT | PyBUF_ND;
2053 if (values == Py_None) {
2054 if (self->substitution_matrix.obj)
2055 PyBuffer_Release(&self->substitution_matrix);
2056 return 0;
2057 }
2058 if (PyObject_GetBuffer(values, &view, flag) != 0) {
2059 PyErr_SetString(PyExc_ValueError, "expected a matrix");
2060 return -1;
2061 }
2062 if (view.ndim != 2) {
2063 PyErr_Format(PyExc_ValueError,
2064 "substitution matrix has incorrect rank (%d expected 2)",
2065 view.ndim);
2066 PyBuffer_Release(&view);
2067 return -1;
2068 }
2069 if (view.len == 0) {
2070 PyErr_SetString(PyExc_ValueError, "substitution matrix has zero size");
2071 PyBuffer_Release(&view);
2072 return -1;
2073 }
2074 if (strcmp(view.format, "d") != 0) {
2075 PyErr_SetString(PyExc_ValueError,
2076 "substitution matrix should contain float values");
2077 PyBuffer_Release(&view);
2078 return -1;
2079 }
2080 if (view.itemsize != sizeof(double)) {
2081 PyErr_Format(PyExc_RuntimeError,
2082 "substitution matrix has unexpected item byte size "
2083 "(%zd, expected %zd)", view.itemsize, sizeof(double));
2084 PyBuffer_Release(&view);
2085 return -1;
2086 }
2087 if (view.shape[0] != view.shape[1]) {
2088 PyErr_Format(PyExc_ValueError,
2089 "substitution matrix should be square "
2090 "(found a %zd x %zd matrix)",
2091 view.shape[0], view.shape[1]);
2092 PyBuffer_Release(&view);
2093 return -1;
2094 }
2095 alphabet = PyObject_GetAttrString(values, "alphabet");
2096 if (alphabet) {
2097 size = set_alphabet(self, alphabet);
2098 Py_DECREF(alphabet);
2099 } else {
2100 /* Set a substitution matrix without setting an alphabet; useful
2101 * when aligning integers. */
2102 PyErr_Clear();
2103 size = set_alphabet(self, Py_None);
2104 }
2105 if (size < 0) {
2106 PyBuffer_Release(&view);
2107 return -1;
2108 }
2109 if (self->substitution_matrix.obj) PyBuffer_Release(&self->substitution_matrix);
2110 self->substitution_matrix = view;
2111 return 0;
2112 }
2113
2114 static char Aligner_alphabet__doc__[] = "alphabet";
2115
2116 static PyObject*
Aligner_get_alphabet(Aligner * self,void * closure)2117 Aligner_get_alphabet(Aligner* self, void* closure)
2118 { PyObject* object = self->alphabet;
2119 if (!object) object = Py_None;
2120 Py_INCREF(object);
2121 return object;
2122 }
2123
2124 static int
Aligner_set_alphabet(Aligner * self,PyObject * alphabet,void * closure)2125 Aligner_set_alphabet(Aligner* self, PyObject* alphabet, void* closure)
2126 {
2127 if (self->substitution_matrix.obj) {
2128 PyErr_SetString(PyExc_AttributeError,
2129 "can't set alphabet if a substitution matrix is used");
2130 return -1;
2131 }
2132 if (set_alphabet(self, alphabet) < 0) return -1;
2133 return 0;
2134 }
2135
2136 static char Aligner_gap_score__doc__[] = "gap score";
2137
2138 static PyObject*
Aligner_get_gap_score(Aligner * self,void * closure)2139 Aligner_get_gap_score(Aligner* self, void* closure)
2140 {
2141 if (self->target_gap_function || self->query_gap_function) {
2142 if (self->target_gap_function != self->query_gap_function) {
2143 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2144 return NULL;
2145 }
2146 Py_INCREF(self->target_gap_function);
2147 return self->target_gap_function;
2148 }
2149 else {
2150 const double score = self->target_internal_open_gap_score;
2151 if (score != self->target_internal_extend_gap_score
2152 || score != self->target_left_open_gap_score
2153 || score != self->target_left_extend_gap_score
2154 || score != self->target_right_open_gap_score
2155 || score != self->target_right_extend_gap_score
2156 || score != self->query_internal_open_gap_score
2157 || score != self->query_internal_extend_gap_score
2158 || score != self->query_left_open_gap_score
2159 || score != self->query_left_extend_gap_score
2160 || score != self->query_right_open_gap_score
2161 || score != self->query_right_extend_gap_score) {
2162 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2163 return NULL;
2164 }
2165 return PyFloat_FromDouble(score);
2166 }
2167 }
2168
2169 static int
Aligner_set_gap_score(Aligner * self,PyObject * value,void * closure)2170 Aligner_set_gap_score(Aligner* self, PyObject* value, void* closure)
2171 { if (PyCallable_Check(value)) {
2172 Py_XDECREF(self->target_gap_function);
2173 Py_XDECREF(self->query_gap_function);
2174 Py_INCREF(value);
2175 Py_INCREF(value);
2176 self->target_gap_function = value;
2177 self->query_gap_function = value;
2178 }
2179 else {
2180 const double score = PyFloat_AsDouble(value);
2181 if (PyErr_Occurred()) return -1;
2182 if (self->target_gap_function) {
2183 Py_DECREF(self->target_gap_function);
2184 self->target_gap_function = NULL;
2185 }
2186 if (self->query_gap_function) {
2187 Py_DECREF(self->query_gap_function);
2188 self->query_gap_function = NULL;
2189 }
2190 self->target_internal_open_gap_score = score;
2191 self->target_internal_extend_gap_score = score;
2192 self->target_left_open_gap_score = score;
2193 self->target_left_extend_gap_score = score;
2194 self->target_right_open_gap_score = score;
2195 self->target_right_extend_gap_score = score;
2196 self->query_internal_open_gap_score = score;
2197 self->query_internal_extend_gap_score = score;
2198 self->query_left_open_gap_score = score;
2199 self->query_left_extend_gap_score = score;
2200 self->query_right_open_gap_score = score;
2201 self->query_right_extend_gap_score = score;
2202 }
2203 self->algorithm = Unknown;
2204 return 0;
2205 }
2206
2207 static char Aligner_open_gap_score__doc__[] = "internal and end open gap score";
2208
2209 static PyObject*
Aligner_get_open_gap_score(Aligner * self,void * closure)2210 Aligner_get_open_gap_score(Aligner* self, void* closure)
2211 {
2212 if (self->target_gap_function || self->query_gap_function) {
2213 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2214 return NULL;
2215 }
2216 else {
2217 const double score = self->target_internal_open_gap_score;
2218 if (score != self->target_left_open_gap_score
2219 || score != self->target_right_open_gap_score
2220 || score != self->query_internal_open_gap_score
2221 || score != self->query_left_open_gap_score
2222 || score != self->query_right_open_gap_score) {
2223 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2224 return NULL;
2225 }
2226 return PyFloat_FromDouble(score);
2227 }
2228 }
2229
2230 static int
Aligner_set_open_gap_score(Aligner * self,PyObject * value,void * closure)2231 Aligner_set_open_gap_score(Aligner* self, PyObject* value, void* closure)
2232 { const double score = PyFloat_AsDouble(value);
2233 if (PyErr_Occurred()) return -1;
2234 if (self->target_gap_function) {
2235 Py_DECREF(self->target_gap_function);
2236 self->target_gap_function = NULL;
2237 }
2238 if (self->query_gap_function) {
2239 Py_DECREF(self->query_gap_function);
2240 self->query_gap_function = NULL;
2241 }
2242 self->target_internal_open_gap_score = score;
2243 self->target_left_open_gap_score = score;
2244 self->target_right_open_gap_score = score;
2245 self->query_internal_open_gap_score = score;
2246 self->query_left_open_gap_score = score;
2247 self->query_right_open_gap_score = score;
2248 self->algorithm = Unknown;
2249 return 0;
2250 }
2251
2252 static char Aligner_extend_gap_score__doc__[] = "extend gap score";
2253
2254 static PyObject*
Aligner_get_extend_gap_score(Aligner * self,void * closure)2255 Aligner_get_extend_gap_score(Aligner* self, void* closure)
2256 {
2257 if (self->target_gap_function || self->query_gap_function) {
2258 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2259 return NULL;
2260 }
2261 else {
2262 const double score = self->target_internal_extend_gap_score;
2263 if (score != self->target_left_extend_gap_score
2264 || score != self->target_right_extend_gap_score
2265 || score != self->query_internal_extend_gap_score
2266 || score != self->query_left_extend_gap_score
2267 || score != self->query_right_extend_gap_score) {
2268 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2269 return NULL;
2270 }
2271 return PyFloat_FromDouble(score);
2272 }
2273 }
2274
2275 static int
Aligner_set_extend_gap_score(Aligner * self,PyObject * value,void * closure)2276 Aligner_set_extend_gap_score(Aligner* self, PyObject* value, void* closure)
2277 { const double score = PyFloat_AsDouble(value);
2278 if (PyErr_Occurred()) return -1;
2279 if (self->target_gap_function) {
2280 Py_DECREF(self->target_gap_function);
2281 self->target_gap_function = NULL;
2282 }
2283 if (self->query_gap_function) {
2284 Py_DECREF(self->query_gap_function);
2285 self->query_gap_function = NULL;
2286 }
2287 self->target_internal_extend_gap_score = score;
2288 self->target_left_extend_gap_score = score;
2289 self->target_right_extend_gap_score = score;
2290 self->query_internal_extend_gap_score = score;
2291 self->query_left_extend_gap_score = score;
2292 self->query_right_extend_gap_score = score;
2293 self->algorithm = Unknown;
2294 return 0;
2295 }
2296
2297 static char Aligner_internal_gap_score__doc__[] = "internal gap score";
2298
2299 static PyObject*
Aligner_get_internal_gap_score(Aligner * self,void * closure)2300 Aligner_get_internal_gap_score(Aligner* self, void* closure)
2301 { if (self->target_gap_function || self->query_gap_function) {
2302 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2303 return NULL;
2304 }
2305 else {
2306 const double score = self->target_internal_open_gap_score;
2307 if (score != self->target_internal_extend_gap_score
2308 || score != self->query_internal_open_gap_score
2309 || score != self->query_internal_extend_gap_score) {
2310 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2311 return NULL;
2312 }
2313 return PyFloat_FromDouble(score);
2314 }
2315 }
2316
2317 static int
Aligner_set_internal_gap_score(Aligner * self,PyObject * value,void * closure)2318 Aligner_set_internal_gap_score(Aligner* self, PyObject* value, void* closure)
2319 { const double score = PyFloat_AsDouble(value);
2320 if (PyErr_Occurred()) return -1;
2321 if (self->target_gap_function) {
2322 Py_DECREF(self->target_gap_function);
2323 self->target_gap_function = NULL;
2324 }
2325 if (self->query_gap_function) {
2326 Py_DECREF(self->query_gap_function);
2327 self->query_gap_function = NULL;
2328 }
2329 self->target_internal_open_gap_score = score;
2330 self->target_internal_extend_gap_score = score;
2331 self->query_internal_open_gap_score = score;
2332 self->query_internal_extend_gap_score = score;
2333 self->algorithm = Unknown;
2334 return 0;
2335 }
2336
2337 static char Aligner_internal_open_gap_score__doc__[] = "internal open gap score";
2338
2339 static PyObject*
Aligner_get_internal_open_gap_score(Aligner * self,void * closure)2340 Aligner_get_internal_open_gap_score(Aligner* self, void* closure)
2341 { if (self->target_gap_function || self->query_gap_function) {
2342 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2343 return NULL;
2344 }
2345 else {
2346 const double score = self->target_internal_open_gap_score;
2347 if (score != self->query_internal_open_gap_score) {
2348 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2349 return NULL;
2350 }
2351 return PyFloat_FromDouble(score);
2352 }
2353 }
2354
2355 static int
Aligner_set_internal_open_gap_score(Aligner * self,PyObject * value,void * closure)2356 Aligner_set_internal_open_gap_score(Aligner* self, PyObject* value, void* closure)
2357 { const double score = PyFloat_AsDouble(value);
2358 if (PyErr_Occurred()) return -1;
2359 if (self->target_gap_function) {
2360 Py_DECREF(self->target_gap_function);
2361 self->target_gap_function = NULL;
2362 }
2363 if (self->query_gap_function) {
2364 Py_DECREF(self->query_gap_function);
2365 self->query_gap_function = NULL;
2366 }
2367 self->target_internal_open_gap_score = score;
2368 self->query_internal_open_gap_score = score;
2369 self->algorithm = Unknown;
2370 return 0;
2371 }
2372
2373 static char Aligner_internal_extend_gap_score__doc__[] = "internal extend gap score";
2374
2375 static PyObject*
Aligner_get_internal_extend_gap_score(Aligner * self,void * closure)2376 Aligner_get_internal_extend_gap_score(Aligner* self, void* closure)
2377 { if (self->target_gap_function || self->query_gap_function) {
2378 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2379 return NULL;
2380 }
2381 else {
2382 const double score = self->target_internal_extend_gap_score;
2383 if (score != self->query_internal_extend_gap_score) {
2384 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2385 return NULL;
2386 }
2387 return PyFloat_FromDouble(score);
2388 }
2389 }
2390
2391 static int
Aligner_set_internal_extend_gap_score(Aligner * self,PyObject * value,void * closure)2392 Aligner_set_internal_extend_gap_score(Aligner* self, PyObject* value,
2393 void* closure)
2394 { const double score = PyFloat_AsDouble(value);
2395 if (PyErr_Occurred()) return -1;
2396 if (self->target_gap_function) {
2397 Py_DECREF(self->target_gap_function);
2398 self->target_gap_function = NULL;
2399 }
2400 if (self->query_gap_function) {
2401 Py_DECREF(self->query_gap_function);
2402 self->query_gap_function = NULL;
2403 }
2404 self->target_internal_extend_gap_score = score;
2405 self->query_internal_extend_gap_score = score;
2406 self->algorithm = Unknown;
2407 return 0;
2408 }
2409
2410 static char Aligner_end_gap_score__doc__[] = "end gap score";
2411
2412 static PyObject*
Aligner_get_end_gap_score(Aligner * self,void * closure)2413 Aligner_get_end_gap_score(Aligner* self, void* closure)
2414 { if (self->target_gap_function || self->query_gap_function) {
2415 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2416 return NULL;
2417 }
2418 else {
2419 const double score = self->target_left_open_gap_score;
2420 if (score != self->target_left_extend_gap_score
2421 || score != self->target_right_open_gap_score
2422 || score != self->target_right_extend_gap_score
2423 || score != self->query_left_open_gap_score
2424 || score != self->query_left_extend_gap_score
2425 || score != self->query_right_open_gap_score
2426 || score != self->query_right_extend_gap_score) {
2427 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2428 return NULL;
2429 }
2430 return PyFloat_FromDouble(score);
2431 }
2432 }
2433
2434 static int
Aligner_set_end_gap_score(Aligner * self,PyObject * value,void * closure)2435 Aligner_set_end_gap_score(Aligner* self, PyObject* value, void* closure)
2436 { const double score = PyFloat_AsDouble(value);
2437 if (PyErr_Occurred()) return -1;
2438 if (self->target_gap_function) {
2439 Py_DECREF(self->target_gap_function);
2440 self->target_gap_function = NULL;
2441 }
2442 if (self->query_gap_function) {
2443 Py_DECREF(self->query_gap_function);
2444 self->query_gap_function = NULL;
2445 }
2446 self->target_left_open_gap_score = score;
2447 self->target_left_extend_gap_score = score;
2448 self->target_right_open_gap_score = score;
2449 self->target_right_extend_gap_score = score;
2450 self->query_left_open_gap_score = score;
2451 self->query_left_extend_gap_score = score;
2452 self->query_right_open_gap_score = score;
2453 self->query_right_extend_gap_score = score;
2454 self->algorithm = Unknown;
2455 return 0;
2456 }
2457
2458 static char Aligner_end_open_gap_score__doc__[] = "end open gap score";
2459
2460 static PyObject*
Aligner_get_end_open_gap_score(Aligner * self,void * closure)2461 Aligner_get_end_open_gap_score(Aligner* self, void* closure)
2462 { if (self->target_gap_function || self->query_gap_function) {
2463 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2464 return NULL;
2465 }
2466 else {
2467 const double score = self->target_left_open_gap_score;
2468 if (score != self->target_right_open_gap_score
2469 || score != self->query_left_open_gap_score
2470 || score != self->query_right_open_gap_score) {
2471 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2472 return NULL;
2473 }
2474 return PyFloat_FromDouble(score);
2475 }
2476 }
2477
2478 static int
Aligner_set_end_open_gap_score(Aligner * self,PyObject * value,void * closure)2479 Aligner_set_end_open_gap_score(Aligner* self, PyObject* value, void* closure)
2480 { const double score = PyFloat_AsDouble(value);
2481 if (PyErr_Occurred()) return -1;
2482 if (self->target_gap_function) {
2483 Py_DECREF(self->target_gap_function);
2484 self->target_gap_function = NULL;
2485 }
2486 if (self->query_gap_function) {
2487 Py_DECREF(self->query_gap_function);
2488 self->query_gap_function = NULL;
2489 }
2490 self->target_left_open_gap_score = score;
2491 self->target_right_open_gap_score = score;
2492 self->query_left_open_gap_score = score;
2493 self->query_right_open_gap_score = score;
2494 self->algorithm = Unknown;
2495 return 0;
2496 }
2497
2498 static char Aligner_end_extend_gap_score__doc__[] = "end extend gap score";
2499
2500 static PyObject*
Aligner_get_end_extend_gap_score(Aligner * self,void * closure)2501 Aligner_get_end_extend_gap_score(Aligner* self, void* closure)
2502 { if (self->target_gap_function || self->query_gap_function) {
2503 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2504 return NULL;
2505 }
2506 else {
2507 const double score = self->target_left_extend_gap_score;
2508 if (score != self->target_right_extend_gap_score
2509 || score != self->query_left_extend_gap_score
2510 || score != self->query_right_extend_gap_score) {
2511 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2512 return NULL;
2513 }
2514 return PyFloat_FromDouble(score);
2515 }
2516 }
2517
2518 static int
Aligner_set_end_extend_gap_score(Aligner * self,PyObject * value,void * closure)2519 Aligner_set_end_extend_gap_score(Aligner* self, PyObject* value, void* closure)
2520 { const double score = PyFloat_AsDouble(value);
2521 if (PyErr_Occurred()) return -1;
2522 if (self->target_gap_function) {
2523 Py_DECREF(self->target_gap_function);
2524 self->target_gap_function = NULL;
2525 }
2526 if (self->query_gap_function) {
2527 Py_DECREF(self->query_gap_function);
2528 self->query_gap_function = NULL;
2529 }
2530 self->target_left_extend_gap_score = score;
2531 self->target_right_extend_gap_score = score;
2532 self->query_left_extend_gap_score = score;
2533 self->query_right_extend_gap_score = score;
2534 self->algorithm = Unknown;
2535 return 0;
2536 }
2537
2538 static char Aligner_left_gap_score__doc__[] = "left gap score";
2539
2540 static PyObject*
Aligner_get_left_gap_score(Aligner * self,void * closure)2541 Aligner_get_left_gap_score(Aligner* self, void* closure)
2542 { if (self->target_gap_function || self->query_gap_function) {
2543 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2544 return NULL;
2545 }
2546 else {
2547 const double score = self->target_left_open_gap_score;
2548 if (score != self->target_left_extend_gap_score
2549 || score != self->query_left_open_gap_score
2550 || score != self->query_left_extend_gap_score) {
2551 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2552 return NULL;
2553 }
2554 return PyFloat_FromDouble(score);
2555 }
2556 }
2557
2558 static int
Aligner_set_left_gap_score(Aligner * self,PyObject * value,void * closure)2559 Aligner_set_left_gap_score(Aligner* self, PyObject* value, void* closure)
2560 { const double score = PyFloat_AsDouble(value);
2561 if (PyErr_Occurred()) return -1;
2562 if (self->target_gap_function) {
2563 Py_DECREF(self->target_gap_function);
2564 self->target_gap_function = NULL;
2565 }
2566 if (self->query_gap_function) {
2567 Py_DECREF(self->query_gap_function);
2568 self->query_gap_function = NULL;
2569 }
2570 self->target_left_open_gap_score = score;
2571 self->target_left_extend_gap_score = score;
2572 self->query_left_open_gap_score = score;
2573 self->query_left_extend_gap_score = score;
2574 self->algorithm = Unknown;
2575 return 0;
2576 }
2577
2578 static char Aligner_right_gap_score__doc__[] = "right gap score";
2579
2580 static PyObject*
Aligner_get_right_gap_score(Aligner * self,void * closure)2581 Aligner_get_right_gap_score(Aligner* self, void* closure)
2582 { if (self->target_gap_function || self->query_gap_function) {
2583 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2584 return NULL;
2585 }
2586 else {
2587 const double score = self->target_right_open_gap_score;
2588 if (score != self->target_right_extend_gap_score
2589 || score != self->query_right_open_gap_score
2590 || score != self->query_right_extend_gap_score) {
2591 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2592 return NULL;
2593 }
2594 return PyFloat_FromDouble(score);
2595 }
2596 }
2597
2598 static int
Aligner_set_right_gap_score(Aligner * self,PyObject * value,void * closure)2599 Aligner_set_right_gap_score(Aligner* self, PyObject* value, void* closure)
2600 { const double score = PyFloat_AsDouble(value);
2601 if (PyErr_Occurred()) return -1;
2602 if (self->target_gap_function) {
2603 Py_DECREF(self->target_gap_function);
2604 self->target_gap_function = NULL;
2605 }
2606 if (self->query_gap_function) {
2607 Py_DECREF(self->query_gap_function);
2608 self->query_gap_function = NULL;
2609 }
2610 self->target_right_open_gap_score = score;
2611 self->target_right_extend_gap_score = score;
2612 self->query_right_open_gap_score = score;
2613 self->query_right_extend_gap_score = score;
2614 self->algorithm = Unknown;
2615 return 0;
2616 }
2617
2618 static char Aligner_left_open_gap_score__doc__[] = "left open gap score";
2619
2620 static PyObject*
Aligner_get_left_open_gap_score(Aligner * self,void * closure)2621 Aligner_get_left_open_gap_score(Aligner* self, void* closure)
2622 { if (self->target_gap_function || self->query_gap_function) {
2623 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2624 return NULL;
2625 }
2626 else {
2627 const double score = self->target_left_open_gap_score;
2628 if (score != self->query_left_open_gap_score) {
2629 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2630 return NULL;
2631 }
2632 return PyFloat_FromDouble(score);
2633 }
2634 }
2635
2636 static int
Aligner_set_left_open_gap_score(Aligner * self,PyObject * value,void * closure)2637 Aligner_set_left_open_gap_score(Aligner* self, PyObject* value, void* closure)
2638 { const double score = PyFloat_AsDouble(value);
2639 if (PyErr_Occurred()) return -1;
2640 if (self->target_gap_function) {
2641 Py_DECREF(self->target_gap_function);
2642 self->target_gap_function = NULL;
2643 }
2644 if (self->query_gap_function) {
2645 Py_DECREF(self->query_gap_function);
2646 self->query_gap_function = NULL;
2647 }
2648 self->target_left_open_gap_score = score;
2649 self->query_left_open_gap_score = score;
2650 self->algorithm = Unknown;
2651 return 0;
2652 }
2653
2654 static char Aligner_left_extend_gap_score__doc__[] = "left extend gap score";
2655
2656 static PyObject*
Aligner_get_left_extend_gap_score(Aligner * self,void * closure)2657 Aligner_get_left_extend_gap_score(Aligner* self, void* closure)
2658 { if (self->target_gap_function || self->query_gap_function) {
2659 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2660 return NULL;
2661 }
2662 else {
2663 const double score = self->target_left_extend_gap_score;
2664 if (score != self->query_left_extend_gap_score) {
2665 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2666 return NULL;
2667 }
2668 return PyFloat_FromDouble(score);
2669 }
2670 }
2671
2672 static int
Aligner_set_left_extend_gap_score(Aligner * self,PyObject * value,void * closure)2673 Aligner_set_left_extend_gap_score(Aligner* self, PyObject* value, void* closure)
2674 { const double score = PyFloat_AsDouble(value);
2675 if (PyErr_Occurred()) return -1;
2676 if (self->target_gap_function) {
2677 Py_DECREF(self->target_gap_function);
2678 self->target_gap_function = NULL;
2679 }
2680 if (self->query_gap_function) {
2681 Py_DECREF(self->query_gap_function);
2682 self->query_gap_function = NULL;
2683 }
2684 self->target_left_extend_gap_score = score;
2685 self->query_left_extend_gap_score = score;
2686 self->algorithm = Unknown;
2687 return 0;
2688 }
2689
2690 static char Aligner_right_open_gap_score__doc__[] = "right open gap score";
2691
2692 static PyObject*
Aligner_get_right_open_gap_score(Aligner * self,void * closure)2693 Aligner_get_right_open_gap_score(Aligner* self, void* closure)
2694 { if (self->target_gap_function || self->query_gap_function) {
2695 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2696 return NULL;
2697 }
2698 else {
2699 const double score = self->target_right_open_gap_score;
2700 if (score != self->query_right_open_gap_score) {
2701 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2702 return NULL;
2703 }
2704 return PyFloat_FromDouble(score);
2705 }
2706 }
2707
2708 static int
Aligner_set_right_open_gap_score(Aligner * self,PyObject * value,void * closure)2709 Aligner_set_right_open_gap_score(Aligner* self, PyObject* value, void* closure)
2710 { const double score = PyFloat_AsDouble(value);
2711 if (PyErr_Occurred()) return -1;
2712 if (self->target_gap_function) {
2713 Py_DECREF(self->target_gap_function);
2714 self->target_gap_function = NULL;
2715 }
2716 if (self->query_gap_function) {
2717 Py_DECREF(self->query_gap_function);
2718 self->query_gap_function = NULL;
2719 }
2720 self->target_right_open_gap_score = score;
2721 self->query_right_open_gap_score = score;
2722 self->algorithm = Unknown;
2723 return 0;
2724 }
2725
2726 static char Aligner_right_extend_gap_score__doc__[] = "right extend gap score";
2727
2728 static PyObject*
Aligner_get_right_extend_gap_score(Aligner * self,void * closure)2729 Aligner_get_right_extend_gap_score(Aligner* self, void* closure)
2730 { if (self->target_gap_function || self->query_gap_function) {
2731 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2732 return NULL;
2733 }
2734 else {
2735 const double score = self->target_right_extend_gap_score;
2736 if (score != self->query_right_extend_gap_score) {
2737 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2738 return NULL;
2739 }
2740 return PyFloat_FromDouble(score);
2741 }
2742 }
2743
2744 static int
Aligner_set_right_extend_gap_score(Aligner * self,PyObject * value,void * closure)2745 Aligner_set_right_extend_gap_score(Aligner* self, PyObject* value, void* closure)
2746 { const double score = PyFloat_AsDouble(value);
2747 if (PyErr_Occurred()) return -1;
2748 if (self->target_gap_function) {
2749 Py_DECREF(self->target_gap_function);
2750 self->target_gap_function = NULL;
2751 }
2752 if (self->query_gap_function) {
2753 Py_DECREF(self->query_gap_function);
2754 self->query_gap_function = NULL;
2755 }
2756 self->target_right_extend_gap_score = score;
2757 self->query_right_extend_gap_score = score;
2758 self->algorithm = Unknown;
2759 return 0;
2760 }
2761
2762 static char Aligner_target_open_gap_score__doc__[] = "target open gap score";
2763
2764 static PyObject*
Aligner_get_target_open_gap_score(Aligner * self,void * closure)2765 Aligner_get_target_open_gap_score(Aligner* self, void* closure)
2766 { if (self->target_gap_function) {
2767 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2768 return NULL;
2769 }
2770 else {
2771 const double score = self->target_internal_open_gap_score;
2772 if (score != self->target_left_open_gap_score
2773 || score != self->target_right_open_gap_score) {
2774 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2775 return NULL;
2776 }
2777 return PyFloat_FromDouble(score);
2778 }
2779 }
2780
2781 static int
Aligner_set_target_open_gap_score(Aligner * self,PyObject * value,void * closure)2782 Aligner_set_target_open_gap_score(Aligner* self, PyObject* value, void* closure)
2783 { const double score = PyFloat_AsDouble(value);
2784 if (PyErr_Occurred()) return -1;
2785 self->target_internal_open_gap_score = score;
2786 self->target_left_open_gap_score = score;
2787 self->target_right_open_gap_score = score;
2788 if (self->target_gap_function) {
2789 Py_DECREF(self->target_gap_function);
2790 self->target_gap_function = NULL;
2791 }
2792 self->algorithm = Unknown;
2793 return 0;
2794 }
2795
2796 static char Aligner_target_extend_gap_score__doc__[] = "target extend gap score";
2797
2798 static PyObject*
Aligner_get_target_extend_gap_score(Aligner * self,void * closure)2799 Aligner_get_target_extend_gap_score(Aligner* self, void* closure)
2800 { if (self->target_gap_function) {
2801 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2802 return NULL;
2803 }
2804 else {
2805 const double score = self->target_internal_extend_gap_score;
2806 if (score != self->target_left_extend_gap_score
2807 || score != self->target_right_extend_gap_score) {
2808 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2809 return NULL;
2810 }
2811 return PyFloat_FromDouble(score);
2812 }
2813 }
2814
2815 static int
Aligner_set_target_extend_gap_score(Aligner * self,PyObject * value,void * closure)2816 Aligner_set_target_extend_gap_score(Aligner* self, PyObject* value, void* closure)
2817 { const double score = PyFloat_AsDouble(value);
2818 if (PyErr_Occurred()) return -1;
2819 self->target_internal_extend_gap_score = score;
2820 self->target_left_extend_gap_score = score;
2821 self->target_right_extend_gap_score = score;
2822 if (self->target_gap_function) {
2823 Py_DECREF(self->target_gap_function);
2824 self->target_gap_function = NULL;
2825 }
2826 self->algorithm = Unknown;
2827 return 0;
2828 }
2829
2830 static char Aligner_target_gap_score__doc__[] = "target gap score";
2831
2832 static PyObject*
Aligner_get_target_gap_score(Aligner * self,void * closure)2833 Aligner_get_target_gap_score(Aligner* self, void* closure)
2834 { if (self->target_gap_function) {
2835 Py_INCREF(self->target_gap_function);
2836 return self->target_gap_function;
2837 }
2838 else {
2839 const double score = self->target_internal_open_gap_score;
2840 if (score != self->target_internal_extend_gap_score
2841 || score != self->target_left_open_gap_score
2842 || score != self->target_left_extend_gap_score
2843 || score != self->target_right_open_gap_score
2844 || score != self->target_right_extend_gap_score) {
2845 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2846 return NULL;
2847 }
2848 return PyFloat_FromDouble(score);
2849 }
2850 }
2851
2852 static int
Aligner_set_target_gap_score(Aligner * self,PyObject * value,void * closure)2853 Aligner_set_target_gap_score(Aligner* self, PyObject* value, void* closure)
2854 {
2855 if (PyCallable_Check(value)) {
2856 Py_XDECREF(self->target_gap_function);
2857 Py_INCREF(value);
2858 self->target_gap_function = value;
2859 }
2860 else {
2861 const double score = PyFloat_AsDouble(value);
2862 if (PyErr_Occurred()) {
2863 PyErr_SetString(PyExc_ValueError,
2864 "gap score should be numerical or callable");
2865 return -1;
2866 }
2867 self->target_internal_open_gap_score = score;
2868 self->target_internal_extend_gap_score = score;
2869 self->target_left_open_gap_score = score;
2870 self->target_left_extend_gap_score = score;
2871 self->target_right_open_gap_score = score;
2872 self->target_right_extend_gap_score = score;
2873 if (self->target_gap_function) {
2874 Py_DECREF(self->target_gap_function);
2875 self->target_gap_function = NULL;
2876 }
2877 }
2878 self->algorithm = Unknown;
2879 return 0;
2880 }
2881
2882 static char Aligner_query_open_gap_score__doc__[] = "query gap open score";
2883
2884 static PyObject*
Aligner_get_query_open_gap_score(Aligner * self,void * closure)2885 Aligner_get_query_open_gap_score(Aligner* self, void* closure)
2886 { if (self->query_gap_function) {
2887 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2888 return NULL;
2889 }
2890 else {
2891 const double score = self->query_internal_open_gap_score;
2892 if (score != self->query_left_open_gap_score
2893 || score != self->query_right_open_gap_score) {
2894 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2895 return NULL;
2896 }
2897 return PyFloat_FromDouble(score);
2898 }
2899 }
2900
2901 static int
Aligner_set_query_open_gap_score(Aligner * self,PyObject * value,void * closure)2902 Aligner_set_query_open_gap_score(Aligner* self, PyObject* value, void* closure)
2903 { const double score = PyFloat_AsDouble(value);
2904 if (PyErr_Occurred()) return -1;
2905 self->query_internal_open_gap_score = score;
2906 self->query_left_open_gap_score = score;
2907 self->query_right_open_gap_score = score;
2908 if (self->query_gap_function) {
2909 Py_DECREF(self->query_gap_function);
2910 self->query_gap_function = NULL;
2911 }
2912 self->algorithm = Unknown;
2913 return 0;
2914 }
2915
2916 static char Aligner_query_extend_gap_score__doc__[] = "query gap extend score";
2917
2918 static PyObject*
Aligner_get_query_extend_gap_score(Aligner * self,void * closure)2919 Aligner_get_query_extend_gap_score(Aligner* self, void* closure)
2920 { if (self->query_gap_function) {
2921 PyErr_SetString(PyExc_ValueError, "using a gap score function");
2922 return NULL;
2923 }
2924 else {
2925 const double score = self->query_internal_extend_gap_score;
2926 if (score != self->query_left_extend_gap_score
2927 || score != self->query_right_extend_gap_score) {
2928 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2929 return NULL;
2930 }
2931 return PyFloat_FromDouble(score);
2932 }
2933 }
2934
2935 static int
Aligner_set_query_extend_gap_score(Aligner * self,PyObject * value,void * closure)2936 Aligner_set_query_extend_gap_score(Aligner* self, PyObject* value, void* closure)
2937 { const double score = PyFloat_AsDouble(value);
2938 if (PyErr_Occurred()) return -1;
2939 self->query_internal_extend_gap_score = score;
2940 self->query_left_extend_gap_score = score;
2941 self->query_right_extend_gap_score = score;
2942 if (self->query_gap_function) {
2943 Py_DECREF(self->query_gap_function);
2944 self->query_gap_function = NULL;
2945 }
2946 self->algorithm = Unknown;
2947 return 0;
2948 }
2949
2950 static char Aligner_query_gap_score__doc__[] = "query gap score";
2951
2952 static PyObject*
Aligner_get_query_gap_score(Aligner * self,void * closure)2953 Aligner_get_query_gap_score(Aligner* self, void* closure)
2954 { if (self->query_gap_function) {
2955 Py_INCREF(self->query_gap_function);
2956 return self->query_gap_function;
2957 }
2958 else {
2959 const double score = self->query_internal_open_gap_score;
2960 if (score != self->query_left_open_gap_score
2961 || score != self->query_right_open_gap_score
2962 || score != self->query_internal_extend_gap_score
2963 || score != self->query_left_extend_gap_score
2964 || score != self->query_right_extend_gap_score) {
2965 PyErr_SetString(PyExc_ValueError, "gap scores are different");
2966 return NULL;
2967 }
2968 return PyFloat_FromDouble(score);
2969 }
2970 }
2971
2972 static int
Aligner_set_query_gap_score(Aligner * self,PyObject * value,void * closure)2973 Aligner_set_query_gap_score(Aligner* self, PyObject* value, void* closure)
2974 { if (PyCallable_Check(value)) {
2975 Py_XDECREF(self->query_gap_function);
2976 Py_INCREF(value);
2977 self->query_gap_function = value;
2978 }
2979 else {
2980 const double score = PyFloat_AsDouble(value);
2981 if (PyErr_Occurred()) {
2982 PyErr_SetString(PyExc_ValueError,
2983 "gap score should be numerical or callable");
2984 return -1;
2985 }
2986 self->query_internal_open_gap_score = score;
2987 self->query_internal_extend_gap_score = score;
2988 self->query_left_open_gap_score = score;
2989 self->query_left_extend_gap_score = score;
2990 self->query_right_open_gap_score = score;
2991 self->query_right_extend_gap_score = score;
2992 if (self->query_gap_function) {
2993 Py_DECREF(self->query_gap_function);
2994 self->query_gap_function = NULL;
2995 }
2996 }
2997 self->algorithm = Unknown;
2998 return 0;
2999 }
3000
3001 static char Aligner_target_internal_open_gap_score__doc__[] = "target internal open gap score";
3002
3003 static PyObject*
Aligner_get_target_internal_open_gap_score(Aligner * self,void * closure)3004 Aligner_get_target_internal_open_gap_score(Aligner* self, void* closure)
3005 { if (self->target_gap_function) {
3006 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3007 return NULL;
3008 }
3009 return PyFloat_FromDouble(self->target_internal_open_gap_score);
3010 }
3011
3012 static int
Aligner_set_target_internal_open_gap_score(Aligner * self,PyObject * value,void * closure)3013 Aligner_set_target_internal_open_gap_score(Aligner* self,
3014 PyObject* value, void* closure)
3015 { const double score = PyFloat_AsDouble(value);
3016 if (PyErr_Occurred()) return -1;
3017 self->target_internal_open_gap_score = score;
3018 if (self->target_gap_function) {
3019 Py_DECREF(self->target_gap_function);
3020 self->target_gap_function = NULL;
3021 }
3022 self->algorithm = Unknown;
3023 return 0;
3024 }
3025
3026 static char Aligner_target_internal_extend_gap_score__doc__[] = "target internal extend gap score";
3027
3028 static PyObject*
Aligner_get_target_internal_extend_gap_score(Aligner * self,void * closure)3029 Aligner_get_target_internal_extend_gap_score(Aligner* self, void* closure)
3030 { if (self->target_gap_function) {
3031 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3032 return NULL;
3033 }
3034 return PyFloat_FromDouble(self->target_internal_extend_gap_score);
3035 }
3036
3037 static int
Aligner_set_target_internal_extend_gap_score(Aligner * self,PyObject * value,void * closure)3038 Aligner_set_target_internal_extend_gap_score(Aligner* self,
3039 PyObject* value, void* closure)
3040 { const double score = PyFloat_AsDouble(value);
3041 if (PyErr_Occurred()) return -1;
3042 self->target_internal_extend_gap_score = score;
3043 if (self->target_gap_function) {
3044 Py_DECREF(self->target_gap_function);
3045 self->target_gap_function = NULL;
3046 }
3047 self->algorithm = Unknown;
3048 return 0;
3049 }
3050
3051 static char Aligner_target_internal_gap_score__doc__[] = "target internal gap score";
3052
3053 static PyObject*
Aligner_get_target_internal_gap_score(Aligner * self,void * closure)3054 Aligner_get_target_internal_gap_score(Aligner* self, void* closure)
3055 { if (self->target_gap_function) {
3056 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3057 return NULL;
3058 }
3059 else {
3060 const double score = self->target_internal_open_gap_score;
3061 if (score != self->target_internal_extend_gap_score) {
3062 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3063 return NULL;
3064 }
3065 return PyFloat_FromDouble(score);
3066 }
3067 }
3068
3069 static int
Aligner_set_target_internal_gap_score(Aligner * self,PyObject * value,void * closure)3070 Aligner_set_target_internal_gap_score(Aligner* self, PyObject* value,
3071 void* closure)
3072 { const double score = PyFloat_AsDouble(value);
3073 if (PyErr_Occurred()) return -1;
3074 self->target_internal_open_gap_score = score;
3075 self->target_internal_extend_gap_score = score;
3076 if (self->target_gap_function) {
3077 Py_DECREF(self->target_gap_function);
3078 self->target_gap_function = NULL;
3079 }
3080 self->algorithm = Unknown;
3081 return 0;
3082 }
3083
3084 static char Aligner_target_end_gap_score__doc__[] = "target end gap score";
3085
3086 static PyObject*
Aligner_get_target_end_gap_score(Aligner * self,void * closure)3087 Aligner_get_target_end_gap_score(Aligner* self, void* closure)
3088 { if (self->target_gap_function) {
3089 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3090 return NULL;
3091 }
3092 else {
3093 const double score = self->target_left_open_gap_score;
3094 if (score != self->target_left_extend_gap_score
3095 || score != self->target_right_open_gap_score
3096 || score != self->target_right_extend_gap_score) {
3097 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3098 return NULL;
3099 }
3100 return PyFloat_FromDouble(score);
3101 }
3102 }
3103
3104 static int
Aligner_set_target_end_gap_score(Aligner * self,PyObject * value,void * closure)3105 Aligner_set_target_end_gap_score(Aligner* self, PyObject* value, void* closure) {
3106 const double score = PyFloat_AsDouble(value);
3107 if (PyErr_Occurred()) return -1;
3108 self->target_left_open_gap_score = score;
3109 self->target_left_extend_gap_score = score;
3110 self->target_right_open_gap_score = score;
3111 self->target_right_extend_gap_score = score;
3112 if (self->target_gap_function) {
3113 Py_DECREF(self->target_gap_function);
3114 self->target_gap_function = NULL;
3115 }
3116 self->algorithm = Unknown;
3117 return 0;
3118 }
3119
3120 static char Aligner_target_end_open_gap_score__doc__[] = "target end open gap score";
3121
3122 static PyObject*
Aligner_get_target_end_open_gap_score(Aligner * self,void * closure)3123 Aligner_get_target_end_open_gap_score(Aligner* self, void* closure)
3124 { if (self->target_gap_function) {
3125 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3126 return NULL;
3127 }
3128 else {
3129 const double score = self->target_left_open_gap_score;
3130 if (score != self->target_right_open_gap_score) {
3131 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3132 return NULL;
3133 }
3134 return PyFloat_FromDouble(score);
3135 }
3136 }
3137
3138 static int
Aligner_set_target_end_open_gap_score(Aligner * self,PyObject * value,void * closure)3139 Aligner_set_target_end_open_gap_score(Aligner* self, PyObject* value,
3140 void* closure)
3141 { const double score = PyFloat_AsDouble(value);
3142 if (PyErr_Occurred()) return -1;
3143 self->target_left_open_gap_score = score;
3144 self->target_right_open_gap_score = score;
3145 if (self->target_gap_function) {
3146 Py_DECREF(self->target_gap_function);
3147 self->target_gap_function = NULL;
3148 }
3149 self->algorithm = Unknown;
3150 return 0;
3151 }
3152
3153 static char Aligner_target_end_extend_gap_score__doc__[] = "target end extend gap score";
3154
3155 static PyObject*
Aligner_get_target_end_extend_gap_score(Aligner * self,void * closure)3156 Aligner_get_target_end_extend_gap_score(Aligner* self, void* closure)
3157 { if (self->target_gap_function) {
3158 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3159 return NULL;
3160 }
3161 else {
3162 const double score = self->target_left_extend_gap_score;
3163 if (score != self->target_right_extend_gap_score) {
3164 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3165 return NULL;
3166 }
3167 return PyFloat_FromDouble(score);
3168 }
3169 }
3170
3171 static int
Aligner_set_target_end_extend_gap_score(Aligner * self,PyObject * value,void * closure)3172 Aligner_set_target_end_extend_gap_score(Aligner* self, PyObject* value, void* closure)
3173 { const double score = PyFloat_AsDouble(value);
3174 if (PyErr_Occurred()) return -1;
3175 self->target_left_extend_gap_score = score;
3176 self->target_right_extend_gap_score = score;
3177 if (self->target_gap_function) {
3178 Py_DECREF(self->target_gap_function);
3179 self->target_gap_function = NULL;
3180 }
3181 self->algorithm = Unknown;
3182 return 0;
3183 }
3184
3185 static char Aligner_target_left_open_gap_score__doc__[] = "target left open score";
3186
3187 static PyObject*
Aligner_get_target_left_open_gap_score(Aligner * self,void * closure)3188 Aligner_get_target_left_open_gap_score(Aligner* self, void* closure)
3189 { if (self->target_gap_function) {
3190 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3191 return NULL;
3192 }
3193 return PyFloat_FromDouble(self->target_left_open_gap_score);
3194 }
3195
3196 static int
Aligner_set_target_left_open_gap_score(Aligner * self,PyObject * value,void * closure)3197 Aligner_set_target_left_open_gap_score(Aligner* self, PyObject* value, void* closure)
3198 { const double score = PyFloat_AsDouble(value);
3199 if (PyErr_Occurred()) return -1;
3200 self->target_left_open_gap_score = score;
3201 if (self->target_gap_function) {
3202 Py_DECREF(self->target_gap_function);
3203 self->target_gap_function = NULL;
3204 }
3205 self->algorithm = Unknown;
3206 return 0;
3207 }
3208
3209 static char Aligner_target_left_extend_gap_score__doc__[] = "target left extend score";
3210
3211 static PyObject*
Aligner_get_target_left_extend_gap_score(Aligner * self,void * closure)3212 Aligner_get_target_left_extend_gap_score(Aligner* self, void* closure)
3213 { if (self->target_gap_function) {
3214 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3215 return NULL;
3216 }
3217 return PyFloat_FromDouble(self->target_left_extend_gap_score);
3218 }
3219
3220 static int
Aligner_set_target_left_extend_gap_score(Aligner * self,PyObject * value,void * closure)3221 Aligner_set_target_left_extend_gap_score(Aligner* self, PyObject* value, void* closure)
3222 { const double score = PyFloat_AsDouble(value);
3223 if (PyErr_Occurred()) return -1;
3224 self->target_left_extend_gap_score = score;
3225 if (self->target_gap_function) {
3226 Py_DECREF(self->target_gap_function);
3227 self->target_gap_function = NULL;
3228 }
3229 self->algorithm = Unknown;
3230 return 0;
3231 }
3232
3233 static char Aligner_target_left_gap_score__doc__[] = "target left score";
3234
3235 static PyObject*
Aligner_get_target_left_gap_score(Aligner * self,void * closure)3236 Aligner_get_target_left_gap_score(Aligner* self, void* closure)
3237 { if (self->target_gap_function) {
3238 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3239 return NULL;
3240 }
3241 else {
3242 const double score = self->target_left_open_gap_score;
3243 if (score != self->target_left_extend_gap_score) {
3244 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3245 return NULL;
3246 }
3247 return PyFloat_FromDouble(score);
3248 }
3249 }
3250
3251 static int
Aligner_set_target_left_gap_score(Aligner * self,PyObject * value,void * closure)3252 Aligner_set_target_left_gap_score(Aligner* self, PyObject* value, void* closure)
3253 { const double score = PyFloat_AsDouble(value);
3254 if (PyErr_Occurred()) return -1;
3255 self->target_left_open_gap_score = score;
3256 self->target_left_extend_gap_score = score;
3257 if (self->target_gap_function) {
3258 Py_DECREF(self->target_gap_function);
3259 self->target_gap_function = NULL;
3260 }
3261 self->algorithm = Unknown;
3262 return 0;
3263 }
3264
3265 static char Aligner_target_right_gap_score_open__doc__[] = "target right open score";
3266
3267 static PyObject*
Aligner_get_target_right_open_gap_score(Aligner * self,void * closure)3268 Aligner_get_target_right_open_gap_score(Aligner* self, void* closure)
3269 { if (self->target_gap_function) {
3270 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3271 return NULL;
3272 }
3273 return PyFloat_FromDouble(self->target_right_open_gap_score);
3274 }
3275
3276 static int
Aligner_set_target_right_open_gap_score(Aligner * self,PyObject * value,void * closure)3277 Aligner_set_target_right_open_gap_score(Aligner* self, PyObject* value, void* closure)
3278 { const double score = PyFloat_AsDouble(value);
3279 if (PyErr_Occurred()) return -1;
3280 self->target_right_open_gap_score = score;
3281 if (self->target_gap_function) {
3282 Py_DECREF(self->target_gap_function);
3283 self->target_gap_function = NULL;
3284 }
3285 self->algorithm = Unknown;
3286 return 0;
3287 }
3288
3289 static char Aligner_target_right_extend_gap_score__doc__[] = "target right extend score";
3290
3291 static PyObject*
Aligner_get_target_right_extend_gap_score(Aligner * self,void * closure)3292 Aligner_get_target_right_extend_gap_score(Aligner* self, void* closure)
3293 { if (self->target_gap_function) {
3294 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3295 return NULL;
3296 }
3297 return PyFloat_FromDouble(self->target_right_extend_gap_score);
3298 }
3299
3300 static int
Aligner_set_target_right_extend_gap_score(Aligner * self,PyObject * value,void * closure)3301 Aligner_set_target_right_extend_gap_score(Aligner* self, PyObject* value, void* closure)
3302 { const double score = PyFloat_AsDouble(value);
3303 if (PyErr_Occurred()) return -1;
3304 self->target_right_extend_gap_score = score;
3305 if (self->target_gap_function) {
3306 Py_DECREF(self->target_gap_function);
3307 self->target_gap_function = NULL;
3308 }
3309 self->algorithm = Unknown;
3310 return 0;
3311 }
3312
3313 static char Aligner_target_right_gap_score__doc__[] = "target right score";
3314
3315 static PyObject*
Aligner_get_target_right_gap_score(Aligner * self,void * closure)3316 Aligner_get_target_right_gap_score(Aligner* self, void* closure)
3317 { if (self->target_gap_function) {
3318 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3319 return NULL;
3320 }
3321 else {
3322 const double score = self->target_right_open_gap_score;
3323 if (score != self->target_right_extend_gap_score) {
3324 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3325 return NULL;
3326 }
3327 return PyFloat_FromDouble(score);
3328 }
3329 }
3330
3331 static int
Aligner_set_target_right_gap_score(Aligner * self,PyObject * value,void * closure)3332 Aligner_set_target_right_gap_score(Aligner* self, PyObject* value, void* closure)
3333 { const double score = PyFloat_AsDouble(value);
3334 if (PyErr_Occurred()) return -1;
3335 self->target_right_open_gap_score = score;
3336 self->target_right_extend_gap_score = score;
3337 if (self->target_gap_function) {
3338 Py_DECREF(self->target_gap_function);
3339 self->target_gap_function = NULL;
3340 }
3341 self->algorithm = Unknown;
3342 return 0;
3343 }
3344
3345 static char Aligner_query_end_gap_score__doc__[] = "query end score";
3346
3347 static PyObject*
Aligner_get_query_end_gap_score(Aligner * self,void * closure)3348 Aligner_get_query_end_gap_score(Aligner* self, void* closure)
3349 { if (self->query_gap_function) {
3350 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3351 return NULL;
3352 }
3353 else {
3354 const double score = self->query_left_open_gap_score;
3355 if (score != self->query_left_extend_gap_score
3356 || score != self->query_right_open_gap_score
3357 || score != self->query_right_extend_gap_score) {
3358 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3359 return NULL;
3360 }
3361 return PyFloat_FromDouble(score);
3362 }
3363 }
3364
3365 static int
Aligner_set_query_end_gap_score(Aligner * self,PyObject * value,void * closure)3366 Aligner_set_query_end_gap_score(Aligner* self, PyObject* value, void* closure)
3367 { const double score = PyFloat_AsDouble(value);
3368 if (PyErr_Occurred()) return -1;
3369 self->query_left_open_gap_score = score;
3370 self->query_left_extend_gap_score = score;
3371 self->query_right_open_gap_score = score;
3372 self->query_right_extend_gap_score = score;
3373 if (self->query_gap_function) {
3374 Py_DECREF(self->query_gap_function);
3375 self->query_gap_function = NULL;
3376 }
3377 self->algorithm = Unknown;
3378 return 0;
3379 }
3380
3381 static char Aligner_query_end_open_gap_score__doc__[] = "query end open score";
3382
3383 static PyObject*
Aligner_get_query_end_open_gap_score(Aligner * self,void * closure)3384 Aligner_get_query_end_open_gap_score(Aligner* self, void* closure)
3385 { if (self->query_gap_function) {
3386 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3387 return NULL;
3388 }
3389 else {
3390 const double score = self->query_left_open_gap_score;
3391 if (score != self->query_right_open_gap_score) {
3392 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3393 return NULL;
3394 }
3395 return PyFloat_FromDouble(score);
3396 }
3397 }
3398
3399 static int
Aligner_set_query_end_open_gap_score(Aligner * self,PyObject * value,void * closure)3400 Aligner_set_query_end_open_gap_score(Aligner* self, PyObject* value, void* closure)
3401 { const double score = PyFloat_AsDouble(value);
3402 if (PyErr_Occurred()) return -1;
3403 self->query_left_open_gap_score = score;
3404 self->query_right_open_gap_score = score;
3405 if (self->query_gap_function) {
3406 Py_DECREF(self->query_gap_function);
3407 self->query_gap_function = NULL;
3408 }
3409 self->algorithm = Unknown;
3410 return 0;
3411 }
3412
3413 static char Aligner_query_end_extend_gap_score__doc__[] = "query end extend score";
3414
3415 static PyObject*
Aligner_get_query_end_extend_gap_score(Aligner * self,void * closure)3416 Aligner_get_query_end_extend_gap_score(Aligner* self, void* closure)
3417 { if (self->query_gap_function) {
3418 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3419 return NULL;
3420 }
3421 else {
3422 const double score = self->query_left_extend_gap_score;
3423 if (score != self->query_right_extend_gap_score) {
3424 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3425 return NULL;
3426 }
3427 return PyFloat_FromDouble(score);
3428 }
3429 }
3430
3431 static int
Aligner_set_query_end_extend_gap_score(Aligner * self,PyObject * value,void * closure)3432 Aligner_set_query_end_extend_gap_score(Aligner* self, PyObject* value, void* closure)
3433 { const double score = PyFloat_AsDouble(value);
3434 if (PyErr_Occurred()) return -1;
3435 self->query_left_extend_gap_score = score;
3436 self->query_right_extend_gap_score = score;
3437 if (self->query_gap_function) {
3438 Py_DECREF(self->query_gap_function);
3439 self->query_gap_function = NULL;
3440 }
3441 self->algorithm = Unknown;
3442 return 0;
3443 }
3444
3445 static char Aligner_query_internal_open_gap_score__doc__[] = "query internal open gap score";
3446
3447 static PyObject*
Aligner_get_query_internal_open_gap_score(Aligner * self,void * closure)3448 Aligner_get_query_internal_open_gap_score(Aligner* self, void* closure)
3449 { if (self->query_gap_function) {
3450 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3451 return NULL;
3452 }
3453 return PyFloat_FromDouble(self->query_internal_open_gap_score);
3454 }
3455
3456 static int
Aligner_set_query_internal_open_gap_score(Aligner * self,PyObject * value,void * closure)3457 Aligner_set_query_internal_open_gap_score(Aligner* self, PyObject* value,
3458 void* closure)
3459 { const double score = PyFloat_AsDouble(value);
3460 if (PyErr_Occurred()) return -1;
3461 self->query_internal_open_gap_score = score;
3462 if (self->query_gap_function) {
3463 Py_DECREF(self->query_gap_function);
3464 self->query_gap_function = NULL;
3465 }
3466 self->algorithm = Unknown;
3467 return 0;
3468 }
3469
3470 static char Aligner_query_internal_extend_gap_score__doc__[] = "query internal extend gap score";
3471
3472 static PyObject*
Aligner_get_query_internal_extend_gap_score(Aligner * self,void * closure)3473 Aligner_get_query_internal_extend_gap_score(Aligner* self, void* closure)
3474 { if (self->query_gap_function) {
3475 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3476 return NULL;
3477 }
3478 return PyFloat_FromDouble(self->query_internal_extend_gap_score);
3479 }
3480
3481 static int
Aligner_set_query_internal_extend_gap_score(Aligner * self,PyObject * value,void * closure)3482 Aligner_set_query_internal_extend_gap_score(Aligner* self, PyObject* value,
3483 void* closure)
3484 { const double score = PyFloat_AsDouble(value);
3485 if (PyErr_Occurred()) return -1;
3486 self->query_internal_extend_gap_score = score;
3487 if (self->query_gap_function) {
3488 Py_DECREF(self->query_gap_function);
3489 self->query_gap_function = NULL;
3490 }
3491 self->algorithm = Unknown;
3492 return 0;
3493 }
3494
3495 static char Aligner_query_internal_gap_score__doc__[] = "query internal gap score";
3496
3497 static PyObject*
Aligner_get_query_internal_gap_score(Aligner * self,void * closure)3498 Aligner_get_query_internal_gap_score(Aligner* self, void* closure)
3499 { if (self->query_gap_function) {
3500 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3501 return NULL;
3502 }
3503 else {
3504 const double score = self->query_internal_open_gap_score;
3505 if (score != self->query_internal_extend_gap_score) {
3506 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3507 return NULL;
3508 }
3509 return PyFloat_FromDouble(score);
3510 }
3511 }
3512
3513 static int
Aligner_set_query_internal_gap_score(Aligner * self,PyObject * value,void * closure)3514 Aligner_set_query_internal_gap_score(Aligner* self, PyObject* value,
3515 void* closure)
3516 { const double score = PyFloat_AsDouble(value);
3517 if (PyErr_Occurred()) return -1;
3518 self->query_internal_open_gap_score = score;
3519 self->query_internal_extend_gap_score = score;
3520 if (self->query_gap_function) {
3521 Py_DECREF(self->query_gap_function);
3522 self->query_gap_function = NULL;
3523 }
3524 self->algorithm = Unknown;
3525 return 0;
3526 }
3527
3528 static char Aligner_query_left_open_gap_score__doc__[] = "query left open score";
3529
3530 static PyObject*
Aligner_get_query_left_open_gap_score(Aligner * self,void * closure)3531 Aligner_get_query_left_open_gap_score(Aligner* self, void* closure)
3532 { if (self->query_gap_function) {
3533 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3534 return NULL;
3535 }
3536 return PyFloat_FromDouble(self->query_left_open_gap_score);
3537 }
3538
3539 static int
Aligner_set_query_left_open_gap_score(Aligner * self,PyObject * value,void * closure)3540 Aligner_set_query_left_open_gap_score(Aligner* self, PyObject* value, void* closure)
3541 { const double score = PyFloat_AsDouble(value);
3542 if (PyErr_Occurred()) return -1;
3543 self->query_left_open_gap_score = score;
3544 if (self->query_gap_function) {
3545 Py_DECREF(self->query_gap_function);
3546 self->query_gap_function = NULL;
3547 }
3548 self->algorithm = Unknown;
3549 return 0;
3550 }
3551
3552 static char Aligner_query_left_extend_gap_score__doc__[] = "query left extend score";
3553
3554 static PyObject*
Aligner_get_query_left_extend_gap_score(Aligner * self,void * closure)3555 Aligner_get_query_left_extend_gap_score(Aligner* self, void* closure)
3556 { if (self->query_gap_function) {
3557 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3558 return NULL;
3559 }
3560 return PyFloat_FromDouble(self->query_left_extend_gap_score);
3561 }
3562
3563 static int
Aligner_set_query_left_extend_gap_score(Aligner * self,PyObject * value,void * closure)3564 Aligner_set_query_left_extend_gap_score(Aligner* self, PyObject* value, void* closure)
3565 { const double score = PyFloat_AsDouble(value);
3566 if (PyErr_Occurred()) return -1;
3567 self->query_left_extend_gap_score = score;
3568 if (self->query_gap_function) {
3569 Py_DECREF(self->query_gap_function);
3570 self->query_gap_function = NULL;
3571 }
3572 self->algorithm = Unknown;
3573 return 0;
3574 }
3575
3576 static char Aligner_query_left_gap_score__doc__[] = "query left score";
3577
3578 static PyObject*
Aligner_get_query_left_gap_score(Aligner * self,void * closure)3579 Aligner_get_query_left_gap_score(Aligner* self, void* closure)
3580 { if (self->query_gap_function) {
3581 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3582 return NULL;
3583 }
3584 else {
3585 const double score = self->query_left_open_gap_score;
3586 if (score != self->query_left_extend_gap_score) {
3587 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3588 return NULL;
3589 }
3590 return PyFloat_FromDouble(score);
3591 }
3592 }
3593
3594 static int
Aligner_set_query_left_gap_score(Aligner * self,PyObject * value,void * closure)3595 Aligner_set_query_left_gap_score(Aligner* self, PyObject* value, void* closure)
3596 { const double score = PyFloat_AsDouble(value);
3597 if (PyErr_Occurred()) return -1;
3598 self->query_left_open_gap_score = score;
3599 self->query_left_extend_gap_score = score;
3600 if (self->query_gap_function) {
3601 Py_DECREF(self->query_gap_function);
3602 self->query_gap_function = NULL;
3603 }
3604 self->algorithm = Unknown;
3605 return 0;
3606 }
3607
3608 static char Aligner_query_right_open_gap_score__doc__[] = "query right open score";
3609
3610 static PyObject*
Aligner_get_query_right_open_gap_score(Aligner * self,void * closure)3611 Aligner_get_query_right_open_gap_score(Aligner* self, void* closure)
3612 { if (self->query_gap_function) {
3613 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3614 return NULL;
3615 }
3616 return PyFloat_FromDouble(self->query_right_open_gap_score);
3617 }
3618
3619 static int
Aligner_set_query_right_open_gap_score(Aligner * self,PyObject * value,void * closure)3620 Aligner_set_query_right_open_gap_score(Aligner* self, PyObject* value, void* closure)
3621 { const double score = PyFloat_AsDouble(value);
3622 if (PyErr_Occurred()) return -1;
3623 self->query_right_open_gap_score = score;
3624 if (self->query_gap_function) {
3625 Py_DECREF(self->query_gap_function);
3626 self->query_gap_function = NULL;
3627 }
3628 self->algorithm = Unknown;
3629 return 0;
3630 }
3631
3632 static char Aligner_query_right_extend_gap_score__doc__[] = "query right extend score";
3633
3634 static PyObject*
Aligner_get_query_right_extend_gap_score(Aligner * self,void * closure)3635 Aligner_get_query_right_extend_gap_score(Aligner* self, void* closure)
3636 { if (self->query_gap_function) {
3637 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3638 return NULL;
3639 }
3640 return PyFloat_FromDouble(self->query_right_extend_gap_score);
3641 }
3642
3643 static int
Aligner_set_query_right_extend_gap_score(Aligner * self,PyObject * value,void * closure)3644 Aligner_set_query_right_extend_gap_score(Aligner* self, PyObject* value, void* closure)
3645 { const double score = PyFloat_AsDouble(value);
3646 if (PyErr_Occurred()) return -1;
3647 self->query_right_extend_gap_score = score;
3648 if (self->query_gap_function) {
3649 Py_DECREF(self->query_gap_function);
3650 self->query_gap_function = NULL;
3651 }
3652 self->algorithm = Unknown;
3653 return 0;
3654 }
3655
3656 static char Aligner_query_right_gap_score__doc__[] = "query right score";
3657
3658 static PyObject*
Aligner_get_query_right_gap_score(Aligner * self,void * closure)3659 Aligner_get_query_right_gap_score(Aligner* self, void* closure)
3660 { if (self->query_gap_function) {
3661 PyErr_SetString(PyExc_ValueError, "using a gap score function");
3662 return NULL;
3663 }
3664 else {
3665 const double score = self->query_right_open_gap_score;
3666 if (score != self->query_right_extend_gap_score) {
3667 PyErr_SetString(PyExc_ValueError, "gap scores are different");
3668 return NULL;
3669 }
3670 return PyFloat_FromDouble(score);
3671 }
3672 }
3673
3674 static int
Aligner_set_query_right_gap_score(Aligner * self,PyObject * value,void * closure)3675 Aligner_set_query_right_gap_score(Aligner* self, PyObject* value, void* closure)
3676 { const double score = PyFloat_AsDouble(value);
3677 if (PyErr_Occurred()) return -1;
3678 self->query_right_open_gap_score = score;
3679 self->query_right_extend_gap_score = score;
3680 if (self->query_gap_function) {
3681 Py_DECREF(self->query_gap_function);
3682 self->query_gap_function = NULL;
3683 }
3684 self->algorithm = Unknown;
3685 return 0;
3686 }
3687
3688 static char Aligner_epsilon__doc__[] = "roundoff epsilon";
3689
3690 static PyObject*
Aligner_get_epsilon(Aligner * self,void * closure)3691 Aligner_get_epsilon(Aligner* self, void* closure)
3692 { return PyFloat_FromDouble(self->epsilon);
3693 }
3694
3695 static int
Aligner_set_epsilon(Aligner * self,PyObject * value,void * closure)3696 Aligner_set_epsilon(Aligner* self, PyObject* value, void* closure)
3697 { const double epsilon = PyFloat_AsDouble(value);
3698 if (PyErr_Occurred()) return -1;
3699 self->epsilon = epsilon;
3700 self->algorithm = Unknown;
3701 return 0;
3702 }
3703
3704 static PyObject*
Aligner_get_wildcard(Aligner * self,void * closure)3705 Aligner_get_wildcard(Aligner* self, void* closure)
3706 {
3707 if (self->wildcard == -1) {
3708 Py_INCREF(Py_None);
3709 return Py_None;
3710 }
3711 else {
3712 return PyUnicode_FromKindAndData(PyUnicode_4BYTE_KIND, &self->wildcard, 1);
3713 }
3714 }
3715
3716 static int
Aligner_set_wildcard(Aligner * self,PyObject * value,void * closure)3717 Aligner_set_wildcard(Aligner* self, PyObject* value, void* closure)
3718 {
3719 if (value == Py_None) {
3720 self->wildcard = -1;
3721 return 0;
3722 }
3723 if (!PyUnicode_Check(value)) {
3724 PyErr_SetString(PyExc_TypeError,
3725 "wildcard should be a single character, or None");
3726 return -1;
3727 }
3728 if (PyUnicode_READY(value) == -1) return -1;
3729 if (PyUnicode_GET_LENGTH(value) != 1) {
3730 PyErr_SetString(PyExc_ValueError,
3731 "wildcard should be a single character, or None");
3732 return -1;
3733 }
3734 self->wildcard = PyUnicode_READ_CHAR(value, 0);
3735 return 0;
3736 }
3737
3738 static char Aligner_wildcard__doc__[] = "wildcard character";
3739
_get_algorithm(Aligner * self)3740 static Algorithm _get_algorithm(Aligner* self)
3741 {
3742 Algorithm algorithm = self->algorithm;
3743 if (algorithm == Unknown) {
3744 const double target_gap_open = self->target_internal_open_gap_score;
3745 const double query_gap_open = self->query_internal_open_gap_score;
3746 const double target_gap_extend = self->target_internal_extend_gap_score;
3747 const double query_gap_extend = self->query_internal_extend_gap_score;
3748 const double target_left_open = self->target_left_open_gap_score;
3749 const double target_left_extend = self->target_left_extend_gap_score;
3750 const double query_left_open = self->query_left_open_gap_score;
3751 const double target_right_open = self->target_right_open_gap_score;
3752 const double query_right_open = self->query_right_open_gap_score;
3753 const double target_right_extend = self->target_right_extend_gap_score;
3754 const double query_left_extend = self->query_left_extend_gap_score;
3755 const double query_right_extend = self->query_right_extend_gap_score;
3756 if (self->target_gap_function || self->query_gap_function)
3757 algorithm = WatermanSmithBeyer;
3758 else if (target_gap_open == target_gap_extend
3759 && query_gap_open == query_gap_extend
3760 && target_left_open == target_left_extend
3761 && target_right_open == target_right_extend
3762 && query_left_open == query_left_extend
3763 && query_right_open == query_right_extend)
3764 algorithm = NeedlemanWunschSmithWaterman;
3765 else
3766 algorithm = Gotoh;
3767 self->algorithm = algorithm;
3768 }
3769 return algorithm;
3770 }
3771
3772
3773 static char Aligner_algorithm__doc__[] = "alignment algorithm";
3774
3775 static PyObject*
Aligner_get_algorithm(Aligner * self,void * closure)3776 Aligner_get_algorithm(Aligner* self, void* closure)
3777 {
3778 const char* s = NULL;
3779 const Mode mode = self->mode;
3780 const Algorithm algorithm = _get_algorithm(self);
3781 switch (algorithm) {
3782 case NeedlemanWunschSmithWaterman:
3783 switch (mode) {
3784 case Global:
3785 s = "Needleman-Wunsch";
3786 break;
3787 case Local:
3788 s = "Smith-Waterman";
3789 break;
3790 }
3791 break;
3792 case Gotoh:
3793 switch (mode) {
3794 case Global:
3795 s = "Gotoh global alignment algorithm";
3796 break;
3797 case Local:
3798 s = "Gotoh local alignment algorithm";
3799 break;
3800 }
3801 break;
3802 case WatermanSmithBeyer:
3803 switch (mode) {
3804 case Global:
3805 s = "Waterman-Smith-Beyer global alignment algorithm";
3806 break;
3807 case Local:
3808 s = "Waterman-Smith-Beyer local alignment algorithm";
3809 break;
3810 }
3811 break;
3812 case Unknown:
3813 default:
3814 break;
3815 }
3816 return PyUnicode_FromString(s);
3817 }
3818
3819 static PyGetSetDef Aligner_getset[] = {
3820 {"mode",
3821 (getter)Aligner_get_mode,
3822 (setter)Aligner_set_mode,
3823 Aligner_mode__doc__, NULL},
3824 {"match_score",
3825 (getter)Aligner_get_match_score,
3826 (setter)Aligner_set_match_score,
3827 Aligner_match_score__doc__, NULL},
3828 {"mismatch_score",
3829 (getter)Aligner_get_mismatch_score,
3830 (setter)Aligner_set_mismatch_score,
3831 Aligner_mismatch_score__doc__, NULL},
3832 {"match", /* synonym for match_score */
3833 (getter)Aligner_get_match_score,
3834 (setter)Aligner_set_match_score,
3835 Aligner_match_score__doc__, NULL},
3836 {"mismatch", /* synonym for mismatch_score */
3837 (getter)Aligner_get_mismatch_score,
3838 (setter)Aligner_set_mismatch_score,
3839 Aligner_mismatch_score__doc__, NULL},
3840 {"substitution_matrix",
3841 (getter)Aligner_get_substitution_matrix,
3842 (setter)Aligner_set_substitution_matrix,
3843 Aligner_substitution_matrix__doc__, NULL},
3844 {"alphabet",
3845 (getter)Aligner_get_alphabet,
3846 (setter)Aligner_set_alphabet,
3847 Aligner_alphabet__doc__, NULL},
3848 {"gap_score",
3849 (getter)Aligner_get_gap_score,
3850 (setter)Aligner_set_gap_score,
3851 Aligner_gap_score__doc__, NULL},
3852 {"open_gap_score",
3853 (getter)Aligner_get_open_gap_score,
3854 (setter)Aligner_set_open_gap_score,
3855 Aligner_open_gap_score__doc__, NULL},
3856 {"extend_gap_score",
3857 (getter)Aligner_get_extend_gap_score,
3858 (setter)Aligner_set_extend_gap_score,
3859 Aligner_extend_gap_score__doc__, NULL},
3860 {"internal_gap_score",
3861 (getter)Aligner_get_internal_gap_score,
3862 (setter)Aligner_set_internal_gap_score,
3863 Aligner_internal_gap_score__doc__, NULL},
3864 {"internal_open_gap_score",
3865 (getter)Aligner_get_internal_open_gap_score,
3866 (setter)Aligner_set_internal_open_gap_score,
3867 Aligner_internal_open_gap_score__doc__, NULL},
3868 {"internal_extend_gap_score",
3869 (getter)Aligner_get_internal_extend_gap_score,
3870 (setter)Aligner_set_internal_extend_gap_score,
3871 Aligner_internal_extend_gap_score__doc__, NULL},
3872 {"end_gap_score",
3873 (getter)Aligner_get_end_gap_score,
3874 (setter)Aligner_set_end_gap_score,
3875 Aligner_end_gap_score__doc__, NULL},
3876 {"end_open_gap_score",
3877 (getter)Aligner_get_end_open_gap_score,
3878 (setter)Aligner_set_end_open_gap_score,
3879 Aligner_end_open_gap_score__doc__, NULL},
3880 {"end_extend_gap_score",
3881 (getter)Aligner_get_end_extend_gap_score,
3882 (setter)Aligner_set_end_extend_gap_score,
3883 Aligner_end_extend_gap_score__doc__, NULL},
3884 {"left_gap_score",
3885 (getter)Aligner_get_left_gap_score,
3886 (setter)Aligner_set_left_gap_score,
3887 Aligner_left_gap_score__doc__, NULL},
3888 {"left_open_gap_score",
3889 (getter)Aligner_get_left_open_gap_score,
3890 (setter)Aligner_set_left_open_gap_score,
3891 Aligner_left_open_gap_score__doc__, NULL},
3892 {"left_extend_gap_score",
3893 (getter)Aligner_get_left_extend_gap_score,
3894 (setter)Aligner_set_left_extend_gap_score,
3895 Aligner_left_extend_gap_score__doc__, NULL},
3896 {"right_gap_score",
3897 (getter)Aligner_get_right_gap_score,
3898 (setter)Aligner_set_right_gap_score,
3899 Aligner_right_gap_score__doc__, NULL},
3900 {"right_open_gap_score",
3901 (getter)Aligner_get_right_open_gap_score,
3902 (setter)Aligner_set_right_open_gap_score,
3903 Aligner_right_open_gap_score__doc__, NULL},
3904 {"right_extend_gap_score",
3905 (getter)Aligner_get_right_extend_gap_score,
3906 (setter)Aligner_set_right_extend_gap_score,
3907 Aligner_right_extend_gap_score__doc__, NULL},
3908 {"target_open_gap_score",
3909 (getter)Aligner_get_target_open_gap_score,
3910 (setter)Aligner_set_target_open_gap_score,
3911 Aligner_target_open_gap_score__doc__, NULL},
3912 {"target_extend_gap_score",
3913 (getter)Aligner_get_target_extend_gap_score,
3914 (setter)Aligner_set_target_extend_gap_score,
3915 Aligner_target_extend_gap_score__doc__, NULL},
3916 {"target_gap_score",
3917 (getter)Aligner_get_target_gap_score,
3918 (setter)Aligner_set_target_gap_score,
3919 Aligner_target_gap_score__doc__, NULL},
3920 {"query_open_gap_score",
3921 (getter)Aligner_get_query_open_gap_score,
3922 (setter)Aligner_set_query_open_gap_score,
3923 Aligner_query_open_gap_score__doc__, NULL},
3924 {"query_extend_gap_score",
3925 (getter)Aligner_get_query_extend_gap_score,
3926 (setter)Aligner_set_query_extend_gap_score,
3927 Aligner_query_extend_gap_score__doc__, NULL},
3928 {"query_gap_score",
3929 (getter)Aligner_get_query_gap_score,
3930 (setter)Aligner_set_query_gap_score,
3931 Aligner_query_gap_score__doc__, NULL},
3932 {"target_end_gap_score",
3933 (getter)Aligner_get_target_end_gap_score,
3934 (setter)Aligner_set_target_end_gap_score,
3935 Aligner_target_end_gap_score__doc__, NULL},
3936 {"target_end_open_gap_score",
3937 (getter)Aligner_get_target_end_open_gap_score,
3938 (setter)Aligner_set_target_end_open_gap_score,
3939 Aligner_target_end_open_gap_score__doc__, NULL},
3940 {"target_end_extend_gap_score",
3941 (getter)Aligner_get_target_end_extend_gap_score,
3942 (setter)Aligner_set_target_end_extend_gap_score,
3943 Aligner_target_end_extend_gap_score__doc__, NULL},
3944 {"target_internal_open_gap_score",
3945 (getter)Aligner_get_target_internal_open_gap_score,
3946 (setter)Aligner_set_target_internal_open_gap_score,
3947 Aligner_target_internal_open_gap_score__doc__, NULL},
3948 {"target_internal_extend_gap_score",
3949 (getter)Aligner_get_target_internal_extend_gap_score,
3950 (setter)Aligner_set_target_internal_extend_gap_score,
3951 Aligner_target_internal_extend_gap_score__doc__, NULL},
3952 {"target_internal_gap_score",
3953 (getter)Aligner_get_target_internal_gap_score,
3954 (setter)Aligner_set_target_internal_gap_score,
3955 Aligner_target_internal_gap_score__doc__, NULL},
3956 {"target_left_open_gap_score",
3957 (getter)Aligner_get_target_left_open_gap_score,
3958 (setter)Aligner_set_target_left_open_gap_score,
3959 Aligner_target_left_open_gap_score__doc__, NULL},
3960 {"target_left_extend_gap_score",
3961 (getter)Aligner_get_target_left_extend_gap_score,
3962 (setter)Aligner_set_target_left_extend_gap_score,
3963 Aligner_target_left_extend_gap_score__doc__, NULL},
3964 {"target_left_gap_score",
3965 (getter)Aligner_get_target_left_gap_score,
3966 (setter)Aligner_set_target_left_gap_score,
3967 Aligner_target_left_gap_score__doc__, NULL},
3968 {"target_right_open_gap_score",
3969 (getter)Aligner_get_target_right_open_gap_score,
3970 (setter)Aligner_set_target_right_open_gap_score,
3971 Aligner_target_right_gap_score_open__doc__, NULL},
3972 {"target_right_extend_gap_score",
3973 (getter)Aligner_get_target_right_extend_gap_score,
3974 (setter)Aligner_set_target_right_extend_gap_score,
3975 Aligner_target_right_extend_gap_score__doc__, NULL},
3976 {"target_right_gap_score",
3977 (getter)Aligner_get_target_right_gap_score,
3978 (setter)Aligner_set_target_right_gap_score,
3979 Aligner_target_right_gap_score__doc__, NULL},
3980 {"query_end_gap_score",
3981 (getter)Aligner_get_query_end_gap_score,
3982 (setter)Aligner_set_query_end_gap_score,
3983 Aligner_query_end_gap_score__doc__, NULL},
3984 {"query_end_open_gap_score",
3985 (getter)Aligner_get_query_end_open_gap_score,
3986 (setter)Aligner_set_query_end_open_gap_score,
3987 Aligner_query_end_open_gap_score__doc__, NULL},
3988 {"query_end_extend_gap_score",
3989 (getter)Aligner_get_query_end_extend_gap_score,
3990 (setter)Aligner_set_query_end_extend_gap_score,
3991 Aligner_query_end_extend_gap_score__doc__, NULL},
3992 {"query_internal_open_gap_score",
3993 (getter)Aligner_get_query_internal_open_gap_score,
3994 (setter)Aligner_set_query_internal_open_gap_score,
3995 Aligner_query_internal_open_gap_score__doc__, NULL},
3996 {"query_internal_extend_gap_score",
3997 (getter)Aligner_get_query_internal_extend_gap_score,
3998 (setter)Aligner_set_query_internal_extend_gap_score,
3999 Aligner_query_internal_extend_gap_score__doc__, NULL},
4000 {"query_internal_gap_score",
4001 (getter)Aligner_get_query_internal_gap_score,
4002 (setter)Aligner_set_query_internal_gap_score,
4003 Aligner_query_internal_gap_score__doc__, NULL},
4004 {"query_left_open_gap_score",
4005 (getter)Aligner_get_query_left_open_gap_score,
4006 (setter)Aligner_set_query_left_open_gap_score,
4007 Aligner_query_left_open_gap_score__doc__, NULL},
4008 {"query_left_extend_gap_score",
4009 (getter)Aligner_get_query_left_extend_gap_score,
4010 (setter)Aligner_set_query_left_extend_gap_score,
4011 Aligner_query_left_extend_gap_score__doc__, NULL},
4012 {"query_left_gap_score",
4013 (getter)Aligner_get_query_left_gap_score,
4014 (setter)Aligner_set_query_left_gap_score,
4015 Aligner_query_left_gap_score__doc__, NULL},
4016 {"query_right_open_gap_score",
4017 (getter)Aligner_get_query_right_open_gap_score,
4018 (setter)Aligner_set_query_right_open_gap_score,
4019 Aligner_query_right_open_gap_score__doc__, NULL},
4020 {"query_right_extend_gap_score",
4021 (getter)Aligner_get_query_right_extend_gap_score,
4022 (setter)Aligner_set_query_right_extend_gap_score,
4023 Aligner_query_right_extend_gap_score__doc__, NULL},
4024 {"query_right_gap_score",
4025 (getter)Aligner_get_query_right_gap_score,
4026 (setter)Aligner_set_query_right_gap_score,
4027 Aligner_query_right_gap_score__doc__, NULL},
4028 {"epsilon",
4029 (getter)Aligner_get_epsilon,
4030 (setter)Aligner_set_epsilon,
4031 Aligner_epsilon__doc__, NULL},
4032 {"wildcard",
4033 (getter)Aligner_get_wildcard,
4034 (setter)Aligner_set_wildcard,
4035 Aligner_wildcard__doc__, NULL},
4036 {"algorithm",
4037 (getter)Aligner_get_algorithm,
4038 (setter)NULL,
4039 Aligner_algorithm__doc__, NULL},
4040 {NULL} /* Sentinel */
4041 };
4042
4043 #define SELECT_SCORE_GLOBAL(score1, score2, score3) \
4044 score = score1; \
4045 temp = score2; \
4046 if (temp > score) score = temp; \
4047 temp = score3; \
4048 if (temp > score) score = temp;
4049
4050 #define SELECT_SCORE_WATERMAN_SMITH_BEYER(score1, score2) \
4051 temp = score1 + gapscore; \
4052 if (temp > score) score = temp; \
4053 temp = score2 + gapscore; \
4054 if (temp > score) score = temp;
4055
4056 #define SELECT_SCORE_GOTOH_LOCAL_ALIGN(score1, score2, score3, score4) \
4057 score = score1; \
4058 temp = score2; \
4059 if (temp > score) score = temp; \
4060 temp = score3; \
4061 if (temp > score) score = temp; \
4062 score += score4; \
4063 if (score < 0) score = 0; \
4064 else if (score > maximum) maximum = score;
4065
4066 #define SELECT_SCORE_LOCAL3(score1, score2, score3) \
4067 score = score1; \
4068 temp = score2; \
4069 if (temp > score) score = temp; \
4070 temp = score3; \
4071 if (temp > score) score = temp; \
4072 if (score < 0) score = 0; \
4073 else if (score > maximum) maximum = score;
4074
4075 #define SELECT_SCORE_LOCAL1(score1) \
4076 score = score1; \
4077 if (score < 0) score = 0; \
4078 else if (score > maximum) maximum = score;
4079
4080 #define SELECT_TRACE_NEEDLEMAN_WUNSCH(hgap, vgap, align_score) \
4081 score = temp + (align_score); \
4082 trace = DIAGONAL; \
4083 temp = row[j-1] + hgap; \
4084 if (temp > score + epsilon) { \
4085 score = temp; \
4086 trace = HORIZONTAL; \
4087 } \
4088 else if (temp > score - epsilon) trace |= HORIZONTAL; \
4089 temp = row[j] + vgap; \
4090 if (temp > score + epsilon) { \
4091 score = temp; \
4092 trace = VERTICAL; \
4093 } \
4094 else if (temp > score - epsilon) trace |= VERTICAL; \
4095 temp = row[j]; \
4096 row[j] = score; \
4097 M[i][j].trace = trace;
4098
4099 #define SELECT_TRACE_SMITH_WATERMAN_HVD(align_score) \
4100 trace = DIAGONAL; \
4101 score = temp + (align_score); \
4102 temp = row[j-1] + gap_extend_A; \
4103 if (temp > score + epsilon) { \
4104 score = temp; \
4105 trace = HORIZONTAL; \
4106 } \
4107 else if (temp > score - epsilon) trace |= HORIZONTAL; \
4108 temp = row[j] + gap_extend_B; \
4109 if (temp > score + epsilon) { \
4110 score = temp; \
4111 trace = VERTICAL; \
4112 } \
4113 else if (temp > score - epsilon) trace |= VERTICAL; \
4114 if (score < epsilon) { \
4115 score = 0; \
4116 trace = STARTPOINT; \
4117 } \
4118 else if (trace & DIAGONAL && score > maximum - epsilon) { \
4119 if (score > maximum + epsilon) { \
4120 for ( ; im < i; im++, jm = 0) \
4121 for ( ; jm <= nB; jm++) M[im][jm].trace &= ~ENDPOINT; \
4122 for ( ; jm < j; jm++) M[im][jm].trace &= ~ENDPOINT; \
4123 im = i; \
4124 jm = j; \
4125 } \
4126 trace |= ENDPOINT; \
4127 } \
4128 M[i][j].trace = trace; \
4129 if (score > maximum) maximum = score; \
4130 temp = row[j]; \
4131 row[j] = score;
4132
4133 #define SELECT_TRACE_SMITH_WATERMAN_D(align_score) \
4134 score = temp + (align_score); \
4135 trace = DIAGONAL; \
4136 if (score < epsilon) { \
4137 score = 0; \
4138 } \
4139 else if (trace & DIAGONAL && score > maximum - epsilon) { \
4140 if (score > maximum + epsilon) { \
4141 for ( ; im < i; im++, jm = 0) \
4142 for ( ; jm <= nB; jm++) M[im][jm].trace &= ~ENDPOINT; \
4143 for ( ; jm < j; jm++) M[im][jm].trace &= ~ENDPOINT; \
4144 im = i; \
4145 jm = j; \
4146 } \
4147 trace |= ENDPOINT; \
4148 } \
4149 M[i][j].trace = trace; \
4150 if (score > maximum) maximum = score; \
4151 temp = row[j]; \
4152 row[j] = score
4153
4154 #define SELECT_TRACE_GOTOH_GLOBAL_GAP(matrix, score1, score2, score3) \
4155 trace = M_MATRIX; \
4156 score = score1; \
4157 temp = score2; \
4158 if (temp > score + epsilon) { \
4159 score = temp; \
4160 trace = Ix_MATRIX; \
4161 } \
4162 else if (temp > score - epsilon) trace |= Ix_MATRIX; \
4163 temp = score3; \
4164 if (temp > score + epsilon) { \
4165 score = temp; \
4166 trace = Iy_MATRIX; \
4167 } \
4168 else if (temp > score - epsilon) trace |= Iy_MATRIX; \
4169 gaps[i][j].matrix = trace;
4170
4171 #define SELECT_TRACE_GOTOH_GLOBAL_ALIGN \
4172 trace = M_MATRIX; \
4173 score = M_temp; \
4174 temp = Ix_temp; \
4175 if (temp > score + epsilon) { \
4176 score = Ix_temp; \
4177 trace = Ix_MATRIX; \
4178 } \
4179 else if (temp > score - epsilon) trace |= Ix_MATRIX; \
4180 temp = Iy_temp; \
4181 if (temp > score + epsilon) { \
4182 score = temp; \
4183 trace = Iy_MATRIX; \
4184 } \
4185 else if (temp > score - epsilon) trace |= Iy_MATRIX; \
4186 M[i][j].trace = trace;
4187
4188 #define SELECT_TRACE_GOTOH_LOCAL_ALIGN(align_score) \
4189 trace = M_MATRIX; \
4190 score = M_temp; \
4191 if (Ix_temp > score + epsilon) { \
4192 score = Ix_temp; \
4193 trace = Ix_MATRIX; \
4194 } \
4195 else if (Ix_temp > score - epsilon) trace |= Ix_MATRIX; \
4196 if (Iy_temp > score + epsilon) { \
4197 score = Iy_temp; \
4198 trace = Iy_MATRIX; \
4199 } \
4200 else if (Iy_temp > score - epsilon) trace |= Iy_MATRIX; \
4201 score += (align_score); \
4202 if (score < epsilon) { \
4203 score = 0; \
4204 trace = STARTPOINT; \
4205 } \
4206 else if (score > maximum - epsilon) { \
4207 if (score > maximum + epsilon) { \
4208 maximum = score; \
4209 for ( ; im < i; im++, jm = 0) \
4210 for ( ; jm <= nB; jm++) M[im][jm].trace &= ~ENDPOINT; \
4211 for ( ; jm < j; jm++) M[im][jm].trace &= ~ENDPOINT; \
4212 im = i; \
4213 jm = j; \
4214 } \
4215 trace |= ENDPOINT; \
4216 } \
4217 M[i][j].trace = trace;
4218
4219 #define SELECT_TRACE_GOTOH_LOCAL_GAP(matrix, score1, score2, score3) \
4220 trace = M_MATRIX; \
4221 score = score1; \
4222 temp = score2; \
4223 if (temp > score + epsilon) { \
4224 score = temp; \
4225 trace = Ix_MATRIX; \
4226 } \
4227 else if (temp > score - epsilon) trace |= Ix_MATRIX; \
4228 temp = score3; \
4229 if (temp > score + epsilon) { \
4230 score = temp; \
4231 trace = Iy_MATRIX; \
4232 } \
4233 else if (temp > score - epsilon) trace |= Iy_MATRIX; \
4234 if (score < epsilon) { \
4235 score = -DBL_MAX; \
4236 trace = 0; \
4237 } \
4238 gaps[i][j].matrix = trace;
4239
4240 #define SELECT_TRACE_WATERMAN_SMITH_BEYER_GLOBAL_ALIGN(score4) \
4241 trace = M_MATRIX; \
4242 score = M_row[i-1][j-1]; \
4243 temp = Ix_row[i-1][j-1]; \
4244 if (temp > score + epsilon) { \
4245 score = temp; \
4246 trace = Ix_MATRIX; \
4247 } \
4248 else if (temp > score - epsilon) trace |= Ix_MATRIX; \
4249 temp = Iy_row[i-1][j-1]; \
4250 if (temp > score + epsilon) { \
4251 score = temp; \
4252 trace = Iy_MATRIX; \
4253 } \
4254 else if (temp > score - epsilon) trace |= Iy_MATRIX; \
4255 M_row[i][j] = score + score4; \
4256 M[i][j].trace = trace;
4257
4258 #define SELECT_TRACE_WATERMAN_SMITH_BEYER_GAP(score1, score2) \
4259 temp = score1 + gapscore; \
4260 if (temp > score - epsilon) { \
4261 if (temp > score + epsilon) { \
4262 score = temp; \
4263 nm = 0; \
4264 ng = 0; \
4265 } \
4266 gapM[nm] = gap; \
4267 nm++; \
4268 } \
4269 temp = score2 + gapscore; \
4270 if (temp > score - epsilon) { \
4271 if (temp > score + epsilon) { \
4272 score = temp; \
4273 nm = 0; \
4274 ng = 0; \
4275 } \
4276 gapXY[ng] = gap; \
4277 ng++; \
4278 }
4279
4280 #define SELECT_TRACE_WATERMAN_SMITH_BEYER_ALIGN(score1, score2, score3, score4) \
4281 trace = M_MATRIX; \
4282 score = score1; \
4283 if (score2 > score + epsilon) { \
4284 score = score2; \
4285 trace = Ix_MATRIX; \
4286 } \
4287 else if (score2 > score - epsilon) trace |= Ix_MATRIX; \
4288 if (score3 > score + epsilon) { \
4289 score = score3; \
4290 trace = Iy_MATRIX; \
4291 } \
4292 else if (score3 > score - epsilon) trace |= Iy_MATRIX; \
4293 score += score4; \
4294 if (score < epsilon) { \
4295 score = 0; \
4296 trace = STARTPOINT; \
4297 } \
4298 else if (score > maximum - epsilon) { \
4299 if (score > maximum + epsilon) { \
4300 maximum = score; \
4301 for ( ; im < i; im++, jm = 0) \
4302 for ( ; jm <= nB; jm++) M[im][jm].trace &= ~ENDPOINT; \
4303 for ( ; jm < j; jm++) M[im][jm].trace &= ~ENDPOINT; \
4304 im = i; \
4305 jm = j; \
4306 } \
4307 trace |= ENDPOINT; \
4308 } \
4309 M_row[i][j] = score; \
4310 M[i][j].trace = trace;
4311
4312 /* ----------------- alignment algorithms ----------------- */
4313
4314 #define NEEDLEMANWUNSCH_SCORE(align_score) \
4315 int i; \
4316 int j; \
4317 int kA; \
4318 int kB; \
4319 const double gap_extend_A = self->target_internal_extend_gap_score; \
4320 const double gap_extend_B = self->query_internal_extend_gap_score; \
4321 double score; \
4322 double temp; \
4323 double* row; \
4324 double left_gap_extend_A; \
4325 double right_gap_extend_A; \
4326 double left_gap_extend_B; \
4327 double right_gap_extend_B; \
4328 switch (strand) { \
4329 case '+': \
4330 left_gap_extend_A = self->target_left_extend_gap_score; \
4331 right_gap_extend_A = self->target_right_extend_gap_score; \
4332 left_gap_extend_B = self->query_left_extend_gap_score; \
4333 right_gap_extend_B = self->query_right_extend_gap_score; \
4334 break; \
4335 case '-': \
4336 left_gap_extend_A = self->target_right_extend_gap_score; \
4337 right_gap_extend_A = self->target_left_extend_gap_score; \
4338 left_gap_extend_B = self->query_right_extend_gap_score; \
4339 right_gap_extend_B = self->query_left_extend_gap_score; \
4340 break; \
4341 default: \
4342 PyErr_SetString(PyExc_RuntimeError, "strand was neither '+' nor '-'"); \
4343 return NULL; \
4344 } \
4345 \
4346 /* Needleman-Wunsch algorithm */ \
4347 row = PyMem_Malloc((nB+1)*sizeof(double)); \
4348 if (!row) return PyErr_NoMemory(); \
4349 \
4350 /* The top row of the score matrix is a special case, \
4351 * as there are no previously aligned characters. \
4352 */ \
4353 row[0] = 0.0; \
4354 for (j = 1; j <= nB; j++) row[j] = j * left_gap_extend_A; \
4355 for (i = 1; i < nA; i++) { \
4356 kA = sA[i-1]; \
4357 temp = row[0]; \
4358 row[0] = i * left_gap_extend_B; \
4359 for (j = 1; j < nB; j++) { \
4360 kB = sB[j-1]; \
4361 SELECT_SCORE_GLOBAL(temp + (align_score), \
4362 row[j] + gap_extend_B, \
4363 row[j-1] + gap_extend_A); \
4364 temp = row[j]; \
4365 row[j] = score; \
4366 } \
4367 kB = sB[nB-1]; \
4368 SELECT_SCORE_GLOBAL(temp + (align_score), \
4369 row[nB] + right_gap_extend_B, \
4370 row[nB-1] + gap_extend_A); \
4371 temp = row[nB]; \
4372 row[nB] = score; \
4373 } \
4374 kA = sA[nA-1]; \
4375 temp = row[0]; \
4376 row[0] = nA * right_gap_extend_B; \
4377 for (j = 1; j < nB; j++) { \
4378 kB = sB[j-1]; \
4379 SELECT_SCORE_GLOBAL(temp + (align_score), \
4380 row[j] + gap_extend_B, \
4381 row[j-1] + right_gap_extend_A); \
4382 temp = row[j]; \
4383 row[j] = score; \
4384 } \
4385 kB = sB[nB-1]; \
4386 SELECT_SCORE_GLOBAL(temp + (align_score), \
4387 row[nB] + right_gap_extend_B, \
4388 row[nB-1] + right_gap_extend_A); \
4389 PyMem_Free(row); \
4390 return PyFloat_FromDouble(score);
4391
4392
4393 #define SMITHWATERMAN_SCORE(align_score) \
4394 int i; \
4395 int j; \
4396 int kA; \
4397 int kB; \
4398 const double gap_extend_A = self->target_internal_extend_gap_score; \
4399 const double gap_extend_B = self->query_internal_extend_gap_score; \
4400 double score; \
4401 double* row; \
4402 double temp; \
4403 double maximum = 0; \
4404 \
4405 /* Smith-Waterman algorithm */ \
4406 row = PyMem_Malloc((nB+1)*sizeof(double)); \
4407 if (!row) return PyErr_NoMemory(); \
4408 \
4409 /* The top row of the score matrix is a special case, \
4410 * as there are no previously aligned characters. \
4411 */ \
4412 for (j = 0; j <= nB; j++) \
4413 row[j] = 0; \
4414 for (i = 1; i < nA; i++) { \
4415 kA = sA[i-1]; \
4416 temp = 0; \
4417 for (j = 1; j < nB; j++) { \
4418 kB = sB[j-1]; \
4419 SELECT_SCORE_LOCAL3(temp + (align_score), \
4420 row[j] + gap_extend_B, \
4421 row[j-1] + gap_extend_A); \
4422 temp = row[j]; \
4423 row[j] = score; \
4424 } \
4425 kB = sB[nB-1]; \
4426 SELECT_SCORE_LOCAL1(temp + (align_score)); \
4427 temp = row[nB]; \
4428 row[nB] = score; \
4429 } \
4430 kA = sA[nA-1]; \
4431 temp = 0; \
4432 for (j = 1; j < nB; j++) { \
4433 kB = sB[j-1]; \
4434 SELECT_SCORE_LOCAL1(temp + (align_score)); \
4435 temp = row[j]; \
4436 row[j] = score; \
4437 } \
4438 kB = sB[nB-1]; \
4439 SELECT_SCORE_LOCAL1(temp + (align_score)); \
4440 PyMem_Free(row); \
4441 return PyFloat_FromDouble(maximum);
4442
4443
4444 #define NEEDLEMANWUNSCH_ALIGN(align_score) \
4445 int i; \
4446 int j; \
4447 int kA; \
4448 int kB; \
4449 const double gap_extend_A = self->target_internal_extend_gap_score; \
4450 const double gap_extend_B = self->query_internal_extend_gap_score; \
4451 const double epsilon = self->epsilon; \
4452 Trace** M; \
4453 double score; \
4454 int trace; \
4455 double temp; \
4456 double* row = NULL; \
4457 PathGenerator* paths; \
4458 double left_gap_extend_A; \
4459 double right_gap_extend_A; \
4460 double left_gap_extend_B; \
4461 double right_gap_extend_B; \
4462 switch (strand) { \
4463 case '+': \
4464 left_gap_extend_A = self->target_left_extend_gap_score; \
4465 right_gap_extend_A = self->target_right_extend_gap_score; \
4466 left_gap_extend_B = self->query_left_extend_gap_score; \
4467 right_gap_extend_B = self->query_right_extend_gap_score; \
4468 break; \
4469 case '-': \
4470 left_gap_extend_A = self->target_right_extend_gap_score; \
4471 right_gap_extend_A = self->target_left_extend_gap_score; \
4472 left_gap_extend_B = self->query_right_extend_gap_score; \
4473 right_gap_extend_B = self->query_left_extend_gap_score; \
4474 break; \
4475 default: \
4476 PyErr_SetString(PyExc_RuntimeError, "strand was neither '+' nor '-'"); \
4477 return NULL; \
4478 } \
4479 \
4480 /* Needleman-Wunsch algorithm */ \
4481 paths = PathGenerator_create_NWSW(nA, nB, Global, strand); \
4482 if (!paths) return NULL; \
4483 row = PyMem_Malloc((nB+1)*sizeof(double)); \
4484 if (!row) { \
4485 Py_DECREF(paths); \
4486 return PyErr_NoMemory(); \
4487 } \
4488 M = paths->M; \
4489 row[0] = 0; \
4490 for (j = 1; j <= nB; j++) row[j] = j * left_gap_extend_A; \
4491 for (i = 1; i < nA; i++) { \
4492 temp = row[0]; \
4493 row[0] = i * left_gap_extend_B; \
4494 kA = sA[i-1]; \
4495 for (j = 1; j < nB; j++) { \
4496 kB = sB[j-1]; \
4497 SELECT_TRACE_NEEDLEMAN_WUNSCH(gap_extend_A, gap_extend_B, align_score); \
4498 } \
4499 kB = sB[j-1]; \
4500 SELECT_TRACE_NEEDLEMAN_WUNSCH(gap_extend_A, right_gap_extend_B, align_score); \
4501 } \
4502 temp = row[0]; \
4503 row[0] = i * left_gap_extend_B; \
4504 kA = sA[nA-1]; \
4505 for (j = 1; j < nB; j++) { \
4506 kB = sB[j-1]; \
4507 SELECT_TRACE_NEEDLEMAN_WUNSCH(right_gap_extend_A, gap_extend_B, align_score); \
4508 } \
4509 kB = sB[j-1]; \
4510 SELECT_TRACE_NEEDLEMAN_WUNSCH(right_gap_extend_A, right_gap_extend_B, align_score); \
4511 PyMem_Free(row); \
4512 M[nA][nB].path = 0; \
4513 return Py_BuildValue("fN", score, paths);
4514
4515
4516 #define SMITHWATERMAN_ALIGN(align_score) \
4517 int i; \
4518 int j; \
4519 int im = nA; \
4520 int jm = nB; \
4521 int kA; \
4522 int kB; \
4523 const double gap_extend_A = self->target_internal_extend_gap_score; \
4524 const double gap_extend_B = self->query_internal_extend_gap_score; \
4525 const double epsilon = self->epsilon; \
4526 Trace** M = NULL; \
4527 double maximum = 0; \
4528 double score = 0; \
4529 double* row = NULL; \
4530 double temp; \
4531 int trace; \
4532 PathGenerator* paths = NULL; \
4533 \
4534 /* Smith-Waterman algorithm */ \
4535 paths = PathGenerator_create_NWSW(nA, nB, Local, strand); \
4536 if (!paths) return NULL; \
4537 row = PyMem_Malloc((nB+1)*sizeof(double)); \
4538 if (!row) { \
4539 Py_DECREF(paths); \
4540 return PyErr_NoMemory(); \
4541 } \
4542 M = paths->M; \
4543 for (j = 0; j <= nB; j++) row[j] = 0; \
4544 for (i = 1; i < nA; i++) { \
4545 temp = 0; \
4546 kA = sA[i-1]; \
4547 for (j = 1; j < nB; j++) { \
4548 kB = sB[j-1]; \
4549 SELECT_TRACE_SMITH_WATERMAN_HVD(align_score); \
4550 } \
4551 kB = sB[nB-1]; \
4552 SELECT_TRACE_SMITH_WATERMAN_D(align_score); \
4553 } \
4554 temp = 0; \
4555 kA = sA[nA-1]; \
4556 for (j = 1; j < nB; j++) { \
4557 kB = sB[j-1]; \
4558 SELECT_TRACE_SMITH_WATERMAN_D(align_score); \
4559 } \
4560 kB = sB[nB-1]; \
4561 SELECT_TRACE_SMITH_WATERMAN_D(align_score); \
4562 PyMem_Free(row); \
4563 \
4564 /* As we don't allow zero-score extensions to alignments, \
4565 * we need to remove all traces towards an ENDPOINT. \
4566 * In addition, some points then won't have any path to a STARTPOINT. \
4567 * Here, use path as a temporary variable to indicate if the point \
4568 * is reachable from a STARTPOINT. If it is unreachable, remove all \
4569 * traces from it, and don't allow it to be an ENDPOINT. It may still \
4570 * be a valid STARTPOINT. */ \
4571 for (j = 0; j <= nB; j++) M[0][j].path = 1; \
4572 for (i = 1; i <= nA; i++) { \
4573 M[i][0].path = 1; \
4574 for (j = 1; j <= nB; j++) { \
4575 trace = M[i][j].trace; \
4576 /* Remove traces to unreachable points. */ \
4577 if (!M[i-1][j-1].path) trace &= ~DIAGONAL; \
4578 if (!M[i][j-1].path) trace &= ~HORIZONTAL; \
4579 if (!M[i-1][j].path) trace &= ~VERTICAL; \
4580 if (trace & (STARTPOINT | HORIZONTAL | VERTICAL | DIAGONAL)) { \
4581 /* The point is reachable. */ \
4582 if (trace & ENDPOINT) M[i][j].path = 0; /* no extensions after ENDPOINT */ \
4583 else M[i][j].path = 1; \
4584 } \
4585 else { \
4586 /* The point is not reachable. Then it is not a STARTPOINT, \
4587 * all traces from it can be removed, and it cannot act as \
4588 * an ENDPOINT. */ \
4589 M[i][j].path = 0; \
4590 trace = 0; \
4591 } \
4592 M[i][j].trace = trace; \
4593 } \
4594 } \
4595 if (maximum == 0) M[0][0].path = NONE; \
4596 else M[0][0].path = 0; \
4597 return Py_BuildValue("fN", maximum, paths);
4598
4599
4600 #define GOTOH_GLOBAL_SCORE(align_score) \
4601 int i; \
4602 int j; \
4603 int kA; \
4604 int kB; \
4605 const double gap_open_A = self->target_internal_open_gap_score; \
4606 const double gap_open_B = self->query_internal_open_gap_score; \
4607 const double gap_extend_A = self->target_internal_extend_gap_score; \
4608 const double gap_extend_B = self->query_internal_extend_gap_score; \
4609 double left_gap_open_A; \
4610 double left_gap_open_B; \
4611 double left_gap_extend_A; \
4612 double left_gap_extend_B; \
4613 double right_gap_open_A; \
4614 double right_gap_open_B; \
4615 double right_gap_extend_A; \
4616 double right_gap_extend_B; \
4617 double* M_row = NULL; \
4618 double* Ix_row = NULL; \
4619 double* Iy_row = NULL; \
4620 double score; \
4621 double temp; \
4622 double M_temp; \
4623 double Ix_temp; \
4624 double Iy_temp; \
4625 switch (strand) { \
4626 case '+': \
4627 left_gap_open_A = self->target_left_open_gap_score; \
4628 left_gap_open_B = self->query_left_open_gap_score; \
4629 left_gap_extend_A = self->target_left_extend_gap_score; \
4630 left_gap_extend_B = self->query_left_extend_gap_score; \
4631 right_gap_open_A = self->target_right_open_gap_score; \
4632 right_gap_open_B = self->query_right_open_gap_score; \
4633 right_gap_extend_A = self->target_right_extend_gap_score; \
4634 right_gap_extend_B = self->query_right_extend_gap_score; \
4635 break; \
4636 case '-': \
4637 left_gap_open_A = self->target_right_open_gap_score; \
4638 left_gap_open_B = self->query_right_open_gap_score; \
4639 left_gap_extend_A = self->target_right_extend_gap_score; \
4640 left_gap_extend_B = self->query_right_extend_gap_score; \
4641 right_gap_open_A = self->target_left_open_gap_score; \
4642 right_gap_open_B = self->query_left_open_gap_score; \
4643 right_gap_extend_A = self->target_left_extend_gap_score; \
4644 right_gap_extend_B = self->query_left_extend_gap_score; \
4645 break; \
4646 default: \
4647 PyErr_SetString(PyExc_RuntimeError, "strand was neither '+' nor '-'"); \
4648 return NULL; \
4649 } \
4650 \
4651 /* Gotoh algorithm with three states */ \
4652 M_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4653 if (!M_row) goto exit; \
4654 Ix_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4655 if (!Ix_row) goto exit; \
4656 Iy_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4657 if (!Iy_row) goto exit; \
4658 \
4659 /* The top row of the score matrix is a special case, \
4660 * as there are no previously aligned characters. \
4661 */ \
4662 M_row[0] = 0; \
4663 Ix_row[0] = -DBL_MAX; \
4664 Iy_row[0] = -DBL_MAX; \
4665 for (j = 1; j <= nB; j++) { \
4666 M_row[j] = -DBL_MAX; \
4667 Ix_row[j] = -DBL_MAX; \
4668 Iy_row[j] = left_gap_open_A + left_gap_extend_A * (j-1); \
4669 } \
4670 \
4671 for (i = 1; i < nA; i++) { \
4672 M_temp = M_row[0]; \
4673 Ix_temp = Ix_row[0]; \
4674 Iy_temp = Iy_row[0]; \
4675 M_row[0] = -DBL_MAX; \
4676 Ix_row[0] = left_gap_open_B + left_gap_extend_B * (i-1); \
4677 Iy_row[0] = -DBL_MAX; \
4678 kA = sA[i-1]; \
4679 for (j = 1; j < nB; j++) { \
4680 kB = sB[j-1]; \
4681 SELECT_SCORE_GLOBAL(M_temp, \
4682 Ix_temp, \
4683 Iy_temp); \
4684 M_temp = M_row[j]; \
4685 M_row[j] = score + (align_score); \
4686 SELECT_SCORE_GLOBAL(M_temp + gap_open_B, \
4687 Ix_row[j] + gap_extend_B, \
4688 Iy_row[j] + gap_open_B); \
4689 Ix_temp = Ix_row[j]; \
4690 Ix_row[j] = score; \
4691 SELECT_SCORE_GLOBAL(M_row[j-1] + gap_open_A, \
4692 Ix_row[j-1] + gap_open_A, \
4693 Iy_row[j-1] + gap_extend_A); \
4694 Iy_temp = Iy_row[j]; \
4695 Iy_row[j] = score; \
4696 } \
4697 kB = sB[nB-1]; \
4698 SELECT_SCORE_GLOBAL(M_temp, \
4699 Ix_temp, \
4700 Iy_temp); \
4701 M_temp = M_row[nB]; \
4702 M_row[nB] = score + (align_score); \
4703 SELECT_SCORE_GLOBAL(M_temp + right_gap_open_B, \
4704 Ix_row[nB] + right_gap_extend_B, \
4705 Iy_row[nB] + right_gap_open_B); \
4706 Ix_row[nB] = score; \
4707 SELECT_SCORE_GLOBAL(M_row[nB-1] + gap_open_A, \
4708 Iy_row[nB-1] + gap_extend_A, \
4709 Ix_row[nB-1] + gap_open_A); \
4710 Iy_row[nB] = score; \
4711 } \
4712 \
4713 M_temp = M_row[0]; \
4714 Ix_temp = Ix_row[0]; \
4715 Iy_temp = Iy_row[0]; \
4716 M_row[0] = -DBL_MAX; \
4717 Ix_row[0] = left_gap_open_B + left_gap_extend_B * (i-1); \
4718 Iy_row[0] = -DBL_MAX; \
4719 kA = sA[nA-1]; \
4720 for (j = 1; j < nB; j++) { \
4721 kB = sB[j-1]; \
4722 SELECT_SCORE_GLOBAL(M_temp, \
4723 Ix_temp, \
4724 Iy_temp); \
4725 M_temp = M_row[j]; \
4726 M_row[j] = score + (align_score); \
4727 SELECT_SCORE_GLOBAL(M_temp + gap_open_B, \
4728 Ix_row[j] + gap_extend_B, \
4729 Iy_row[j] + gap_open_B); \
4730 Ix_temp = Ix_row[j]; \
4731 Ix_row[j] = score; \
4732 SELECT_SCORE_GLOBAL(M_row[j-1] + right_gap_open_A, \
4733 Iy_row[j-1] + right_gap_extend_A, \
4734 Ix_row[j-1] + right_gap_open_A); \
4735 Iy_temp = Iy_row[j]; \
4736 Iy_row[j] = score; \
4737 } \
4738 \
4739 kB = sB[nB-1]; \
4740 SELECT_SCORE_GLOBAL(M_temp, \
4741 Ix_temp, \
4742 Iy_temp); \
4743 M_temp = M_row[nB]; \
4744 M_row[nB] = score + (align_score); \
4745 SELECT_SCORE_GLOBAL(M_temp + right_gap_open_B, \
4746 Ix_row[nB] + right_gap_extend_B, \
4747 Iy_row[nB] + right_gap_open_B); \
4748 Ix_temp = Ix_row[nB]; \
4749 Ix_row[nB] = score; \
4750 SELECT_SCORE_GLOBAL(M_row[nB-1] + right_gap_open_A, \
4751 Ix_row[nB-1] + right_gap_open_A, \
4752 Iy_row[nB-1] + right_gap_extend_A); \
4753 Iy_temp = Iy_row[nB]; \
4754 Iy_row[nB] = score; \
4755 \
4756 SELECT_SCORE_GLOBAL(M_row[nB], Ix_row[nB], Iy_row[nB]); \
4757 PyMem_Free(M_row); \
4758 PyMem_Free(Ix_row); \
4759 PyMem_Free(Iy_row); \
4760 return PyFloat_FromDouble(score); \
4761 \
4762 exit: \
4763 if (M_row) PyMem_Free(M_row); \
4764 if (Ix_row) PyMem_Free(Ix_row); \
4765 if (Iy_row) PyMem_Free(Iy_row); \
4766 return PyErr_NoMemory(); \
4767
4768
4769 #define GOTOH_LOCAL_SCORE(align_score) \
4770 int i; \
4771 int j; \
4772 int kA; \
4773 int kB; \
4774 const double gap_open_A = self->target_internal_open_gap_score; \
4775 const double gap_open_B = self->query_internal_open_gap_score; \
4776 const double gap_extend_A = self->target_internal_extend_gap_score; \
4777 const double gap_extend_B = self->query_internal_extend_gap_score; \
4778 double* M_row = NULL; \
4779 double* Ix_row = NULL; \
4780 double* Iy_row = NULL; \
4781 double score; \
4782 double temp; \
4783 double M_temp; \
4784 double Ix_temp; \
4785 double Iy_temp; \
4786 double maximum = 0.0; \
4787 \
4788 /* Gotoh algorithm with three states */ \
4789 M_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4790 if (!M_row) goto exit; \
4791 Ix_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4792 if (!Ix_row) goto exit; \
4793 Iy_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4794 if (!Iy_row) goto exit; \
4795 \
4796 /* The top row of the score matrix is a special case, \
4797 * as there are no previously aligned characters. \
4798 */ \
4799 M_row[0] = 0; \
4800 Ix_row[0] = -DBL_MAX; \
4801 Iy_row[0] = -DBL_MAX; \
4802 for (j = 1; j <= nB; j++) { \
4803 M_row[j] = -DBL_MAX; \
4804 Ix_row[j] = -DBL_MAX; \
4805 Iy_row[j] = 0; \
4806 } \
4807 for (i = 1; i < nA; i++) { \
4808 M_temp = M_row[0]; \
4809 Ix_temp = Ix_row[0]; \
4810 Iy_temp = Iy_row[0]; \
4811 M_row[0] = -DBL_MAX; \
4812 Ix_row[0] = 0; \
4813 Iy_row[0] = -DBL_MAX; \
4814 kA = sA[i-1]; \
4815 for (j = 1; j < nB; j++) { \
4816 kB = sB[j-1]; \
4817 SELECT_SCORE_GOTOH_LOCAL_ALIGN(M_temp, \
4818 Ix_temp, \
4819 Iy_temp, \
4820 (align_score)); \
4821 M_temp = M_row[j]; \
4822 M_row[j] = score; \
4823 SELECT_SCORE_LOCAL3(M_temp + gap_open_B, \
4824 Ix_row[j] + gap_extend_B, \
4825 Iy_row[j] + gap_open_B); \
4826 Ix_temp = Ix_row[j]; \
4827 Ix_row[j] = score; \
4828 SELECT_SCORE_LOCAL3(M_row[j-1] + gap_open_A, \
4829 Ix_row[j-1] + gap_open_A, \
4830 Iy_row[j-1] + gap_extend_A); \
4831 Iy_temp = Iy_row[j]; \
4832 Iy_row[j] = score; \
4833 } \
4834 kB = sB[nB-1]; \
4835 Ix_row[nB] = 0; \
4836 Iy_row[nB] = 0; \
4837 SELECT_SCORE_GOTOH_LOCAL_ALIGN(M_temp, \
4838 Ix_temp, \
4839 Iy_temp, \
4840 (align_score)); \
4841 M_temp = M_row[nB]; \
4842 M_row[nB] = score; \
4843 } \
4844 M_temp = M_row[0]; \
4845 Ix_temp = Ix_row[0]; \
4846 Iy_temp = Iy_row[0]; \
4847 M_row[0] = -DBL_MAX; \
4848 Ix_row[0] = 0; \
4849 Iy_row[0] = -DBL_MAX; \
4850 kA = sA[nA-1]; \
4851 for (j = 1; j < nB; j++) { \
4852 kB = sB[j-1]; \
4853 SELECT_SCORE_GOTOH_LOCAL_ALIGN(M_temp, \
4854 Ix_temp, \
4855 Iy_temp, \
4856 (align_score)); \
4857 M_temp = M_row[j]; \
4858 M_row[j] = score; \
4859 Ix_temp = Ix_row[j]; \
4860 Iy_temp = Iy_row[j]; \
4861 Ix_row[j] = 0; \
4862 Iy_row[j] = 0; \
4863 } \
4864 kB = sB[nB-1]; \
4865 SELECT_SCORE_GOTOH_LOCAL_ALIGN(M_temp, \
4866 Ix_temp, \
4867 Iy_temp, \
4868 (align_score)); \
4869 PyMem_Free(M_row); \
4870 PyMem_Free(Ix_row); \
4871 PyMem_Free(Iy_row); \
4872 return PyFloat_FromDouble(maximum); \
4873 exit: \
4874 if (M_row) PyMem_Free(M_row); \
4875 if (Ix_row) PyMem_Free(Ix_row); \
4876 if (Iy_row) PyMem_Free(Iy_row); \
4877 return PyErr_NoMemory(); \
4878
4879
4880 #define GOTOH_GLOBAL_ALIGN(align_score) \
4881 int i; \
4882 int j; \
4883 int kA; \
4884 int kB; \
4885 const double gap_open_A = self->target_internal_open_gap_score; \
4886 const double gap_open_B = self->query_internal_open_gap_score; \
4887 const double gap_extend_A = self->target_internal_extend_gap_score; \
4888 const double gap_extend_B = self->query_internal_extend_gap_score; \
4889 double left_gap_open_A; \
4890 double left_gap_open_B; \
4891 double left_gap_extend_A; \
4892 double left_gap_extend_B; \
4893 double right_gap_open_A; \
4894 double right_gap_open_B; \
4895 double right_gap_extend_A; \
4896 double right_gap_extend_B; \
4897 const double epsilon = self->epsilon; \
4898 TraceGapsGotoh** gaps = NULL; \
4899 Trace** M = NULL; \
4900 double* M_row = NULL; \
4901 double* Ix_row = NULL; \
4902 double* Iy_row = NULL; \
4903 double score; \
4904 int trace; \
4905 double temp; \
4906 double M_temp; \
4907 double Ix_temp; \
4908 double Iy_temp; \
4909 PathGenerator* paths; \
4910 switch (strand) { \
4911 case '+': \
4912 left_gap_open_A = self->target_left_open_gap_score; \
4913 left_gap_open_B = self->query_left_open_gap_score; \
4914 left_gap_extend_A = self->target_left_extend_gap_score; \
4915 left_gap_extend_B = self->query_left_extend_gap_score; \
4916 right_gap_open_A = self->target_right_open_gap_score; \
4917 right_gap_open_B = self->query_right_open_gap_score; \
4918 right_gap_extend_A = self->target_right_extend_gap_score; \
4919 right_gap_extend_B = self->query_right_extend_gap_score; \
4920 break; \
4921 case '-': \
4922 left_gap_open_A = self->target_right_open_gap_score; \
4923 left_gap_open_B = self->query_right_open_gap_score; \
4924 left_gap_extend_A = self->target_right_extend_gap_score; \
4925 left_gap_extend_B = self->query_right_extend_gap_score; \
4926 right_gap_open_A = self->target_left_open_gap_score; \
4927 right_gap_open_B = self->query_left_open_gap_score; \
4928 right_gap_extend_A = self->target_left_extend_gap_score; \
4929 right_gap_extend_B = self->query_left_extend_gap_score; \
4930 break; \
4931 default: \
4932 PyErr_SetString(PyExc_RuntimeError, "strand was neither '+' nor '-'"); \
4933 return NULL; \
4934 } \
4935 \
4936 /* Gotoh algorithm with three states */ \
4937 paths = PathGenerator_create_Gotoh(nA, nB, Global, strand); \
4938 if (!paths) return NULL; \
4939 M_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4940 if (!M_row) goto exit; \
4941 Ix_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4942 if (!Ix_row) goto exit; \
4943 Iy_row = PyMem_Malloc((nB+1)*sizeof(double)); \
4944 if (!Iy_row) goto exit; \
4945 M = paths->M; \
4946 gaps = paths->gaps.gotoh; \
4947 \
4948 /* Gotoh algorithm with three states */ \
4949 M_row[0] = 0; \
4950 Ix_row[0] = -DBL_MAX; \
4951 Iy_row[0] = -DBL_MAX; \
4952 for (j = 1; j <= nB; j++) { \
4953 M_row[j] = -DBL_MAX; \
4954 Ix_row[j] = -DBL_MAX; \
4955 Iy_row[j] = left_gap_open_A + left_gap_extend_A * (j-1); \
4956 } \
4957 for (i = 1; i < nA; i++) { \
4958 kA = sA[i-1]; \
4959 M_temp = M_row[0]; \
4960 Ix_temp = Ix_row[0]; \
4961 Iy_temp = Iy_row[0]; \
4962 M_row[0] = -DBL_MAX; \
4963 Ix_row[0] = left_gap_open_B + left_gap_extend_B * (i-1); \
4964 Iy_row[0] = -DBL_MAX; \
4965 for (j = 1; j < nB; j++) { \
4966 kB = sB[j-1]; \
4967 SELECT_TRACE_GOTOH_GLOBAL_ALIGN; \
4968 M_temp = M_row[j]; \
4969 M_row[j] = score + (align_score); \
4970 SELECT_TRACE_GOTOH_GLOBAL_GAP(Ix, \
4971 M_temp + gap_open_B, \
4972 Ix_row[j] + gap_extend_B, \
4973 Iy_row[j] + gap_open_B); \
4974 Ix_temp = Ix_row[j]; \
4975 Ix_row[j] = score; \
4976 SELECT_TRACE_GOTOH_GLOBAL_GAP(Iy, \
4977 M_row[j-1] + gap_open_A, \
4978 Ix_row[j-1] + gap_open_A, \
4979 Iy_row[j-1] + gap_extend_A); \
4980 Iy_temp = Iy_row[j]; \
4981 Iy_row[j] = score; \
4982 } \
4983 kB = sB[nB-1]; \
4984 SELECT_TRACE_GOTOH_GLOBAL_ALIGN; \
4985 M_temp = M_row[nB]; \
4986 M_row[nB] = score + (align_score); \
4987 SELECT_TRACE_GOTOH_GLOBAL_GAP(Ix, \
4988 M_temp + right_gap_open_B, \
4989 Ix_row[nB] + right_gap_extend_B, \
4990 Iy_row[nB] + right_gap_open_B); \
4991 Ix_temp = Ix_row[nB]; \
4992 Ix_row[nB] = score; \
4993 SELECT_TRACE_GOTOH_GLOBAL_GAP(Iy, \
4994 M_row[nB-1] + gap_open_A, \
4995 Ix_row[nB-1] + gap_open_A, \
4996 Iy_row[nB-1] + gap_extend_A); \
4997 Iy_temp = Iy_row[nB]; \
4998 Iy_row[nB] = score; \
4999 } \
5000 kA = sA[nA-1]; \
5001 M_temp = M_row[0]; \
5002 Ix_temp = Ix_row[0]; \
5003 Iy_temp = Iy_row[0]; \
5004 M_row[0] = -DBL_MAX; \
5005 Ix_row[0] = left_gap_open_B + left_gap_extend_B * (nA-1); \
5006 Iy_row[0] = -DBL_MAX; \
5007 for (j = 1; j < nB; j++) { \
5008 kB = sB[j-1]; \
5009 SELECT_TRACE_GOTOH_GLOBAL_ALIGN; \
5010 M_temp = M_row[j]; \
5011 M_row[j] = score + (align_score); \
5012 SELECT_TRACE_GOTOH_GLOBAL_GAP(Ix, \
5013 M_temp + gap_open_B, \
5014 Ix_row[j] + gap_extend_B, \
5015 Iy_row[j] + gap_open_B); \
5016 Ix_temp = Ix_row[j]; \
5017 Ix_row[j] = score; \
5018 SELECT_TRACE_GOTOH_GLOBAL_GAP(Iy, \
5019 M_row[j-1] + right_gap_open_A, \
5020 Ix_row[j-1] + right_gap_open_A, \
5021 Iy_row[j-1] + right_gap_extend_A); \
5022 Iy_temp = Iy_row[j]; \
5023 Iy_row[j] = score; \
5024 } \
5025 kB = sB[nB-1]; \
5026 SELECT_TRACE_GOTOH_GLOBAL_ALIGN; \
5027 M_temp = M_row[j]; \
5028 M_row[j] = score + (align_score); \
5029 SELECT_TRACE_GOTOH_GLOBAL_GAP(Ix, \
5030 M_temp + right_gap_open_B, \
5031 Ix_row[j] + right_gap_extend_B, \
5032 Iy_row[j] + right_gap_open_B); \
5033 Ix_row[nB] = score; \
5034 SELECT_TRACE_GOTOH_GLOBAL_GAP(Iy, \
5035 M_row[j-1] + right_gap_open_A, \
5036 Ix_row[j-1] + right_gap_open_A, \
5037 Iy_row[j-1] + right_gap_extend_A); \
5038 Iy_row[nB] = score; \
5039 M[nA][nB].path = 0; \
5040 \
5041 /* traceback */ \
5042 SELECT_SCORE_GLOBAL(M_row[nB], Ix_row[nB], Iy_row[nB]); \
5043 if (M_row[nB] < score - epsilon) M[nA][nB].trace = 0; \
5044 if (Ix_row[nB] < score - epsilon) gaps[nA][nB].Ix = 0; \
5045 if (Iy_row[nB] < score - epsilon) gaps[nA][nB].Iy = 0; \
5046 return Py_BuildValue("fN", score, paths); \
5047 exit: \
5048 Py_DECREF(paths); \
5049 if (M_row) PyMem_Free(M_row); \
5050 if (Ix_row) PyMem_Free(Ix_row); \
5051 if (Iy_row) PyMem_Free(Iy_row); \
5052 return PyErr_NoMemory(); \
5053
5054
5055 #define GOTOH_LOCAL_ALIGN(align_score) \
5056 int i; \
5057 int j; \
5058 int im = nA; \
5059 int jm = nB; \
5060 int kA; \
5061 int kB; \
5062 const double gap_open_A = self->target_internal_open_gap_score; \
5063 const double gap_open_B = self->query_internal_open_gap_score; \
5064 const double gap_extend_A = self->target_internal_extend_gap_score; \
5065 const double gap_extend_B = self->query_internal_extend_gap_score; \
5066 const double epsilon = self->epsilon; \
5067 Trace** M = NULL; \
5068 TraceGapsGotoh** gaps = NULL; \
5069 double* M_row = NULL; \
5070 double* Ix_row = NULL; \
5071 double* Iy_row = NULL; \
5072 double score; \
5073 int trace; \
5074 double temp; \
5075 double M_temp; \
5076 double Ix_temp; \
5077 double Iy_temp; \
5078 double maximum = 0.0; \
5079 PathGenerator* paths; \
5080 \
5081 /* Gotoh algorithm with three states */ \
5082 paths = PathGenerator_create_Gotoh(nA, nB, Local, strand); \
5083 if (!paths) return NULL; \
5084 M = paths->M; \
5085 gaps = paths->gaps.gotoh; \
5086 M_row = PyMem_Malloc((nB+1)*sizeof(double)); \
5087 if (!M_row) goto exit; \
5088 Ix_row = PyMem_Malloc((nB+1)*sizeof(double)); \
5089 if (!Ix_row) goto exit; \
5090 Iy_row = PyMem_Malloc((nB+1)*sizeof(double)); \
5091 if (!Iy_row) goto exit; \
5092 M_row[0] = 0; \
5093 Ix_row[0] = -DBL_MAX; \
5094 Iy_row[0] = -DBL_MAX; \
5095 for (j = 1; j <= nB; j++) { \
5096 M_row[j] = 0; \
5097 Ix_row[j] = -DBL_MAX; \
5098 Iy_row[j] = -DBL_MAX; \
5099 } \
5100 for (i = 1; i < nA; i++) { \
5101 M_temp = M_row[0]; \
5102 Ix_temp = Ix_row[0]; \
5103 Iy_temp = Iy_row[0]; \
5104 M_row[0] = 0; \
5105 Ix_row[0] = -DBL_MAX; \
5106 Iy_row[0] = -DBL_MAX; \
5107 kA = sA[i-1]; \
5108 for (j = 1; j < nB; j++) { \
5109 kB = sB[j-1]; \
5110 SELECT_TRACE_GOTOH_LOCAL_ALIGN(align_score) \
5111 M_temp = M_row[j]; \
5112 M_row[j] = score; \
5113 SELECT_TRACE_GOTOH_LOCAL_GAP(Ix, \
5114 M_temp + gap_open_B, \
5115 Ix_row[j] + gap_extend_B, \
5116 Iy_row[j] + gap_open_B); \
5117 Ix_temp = Ix_row[j]; \
5118 Ix_row[j] = score; \
5119 SELECT_TRACE_GOTOH_LOCAL_GAP(Iy, \
5120 M_row[j-1] + gap_open_A, \
5121 Ix_row[j-1] + gap_open_A, \
5122 Iy_row[j-1] + gap_extend_A); \
5123 Iy_temp = Iy_row[j]; \
5124 Iy_row[j] = score; \
5125 } \
5126 kB = sB[nB-1]; \
5127 SELECT_TRACE_GOTOH_LOCAL_ALIGN(align_score) \
5128 M_temp = M_row[j]; \
5129 M_row[j] = score; \
5130 Ix_temp = Ix_row[nB]; \
5131 Ix_row[nB] = 0; \
5132 gaps[i][nB].Ix = 0; \
5133 Iy_temp = Iy_row[nB]; \
5134 Iy_row[nB] = 0; \
5135 gaps[i][nB].Iy = 0; \
5136 } \
5137 M_temp = M_row[0]; \
5138 M_row[0] = 0; \
5139 M[nA][0].trace = 0; \
5140 Ix_temp = Ix_row[0]; \
5141 Ix_row[0] = -DBL_MAX; \
5142 gaps[nA][0].Ix = 0; \
5143 gaps[nA][0].Iy = 0; \
5144 Iy_temp = Iy_row[0]; \
5145 Iy_row[0] = -DBL_MAX; \
5146 kA = sA[nA-1]; \
5147 for (j = 1; j < nB; j++) { \
5148 kB = sB[j-1]; \
5149 SELECT_TRACE_GOTOH_LOCAL_ALIGN(align_score) \
5150 M_temp = M_row[j]; \
5151 M_row[j] = score; \
5152 Ix_temp = Ix_row[j]; \
5153 Ix_row[j] = 0; \
5154 gaps[nA][j].Ix = 0; \
5155 Iy_temp = Iy_row[j]; \
5156 Iy_row[j] = 0; \
5157 gaps[nA][j].Iy = 0; \
5158 } \
5159 kB = sB[nB-1]; \
5160 SELECT_TRACE_GOTOH_LOCAL_ALIGN(align_score) \
5161 gaps[nA][nB].Ix = 0; \
5162 gaps[nA][nB].Iy = 0; \
5163 \
5164 PyMem_Free(M_row); \
5165 PyMem_Free(Ix_row); \
5166 PyMem_Free(Iy_row); \
5167 \
5168 /* As we don't allow zero-score extensions to alignments, \
5169 * we need to remove all traces towards an ENDPOINT. \
5170 * In addition, some points then won't have any path to a STARTPOINT. \
5171 * Here, use path as a temporary variable to indicate if the point \
5172 * is reachable from a STARTPOINT. If it is unreachable, remove all \
5173 * traces from it, and don't allow it to be an ENDPOINT. It may still \
5174 * be a valid STARTPOINT. */ \
5175 for (j = 0; j <= nB; j++) M[0][j].path = M_MATRIX; \
5176 for (i = 1; i <= nA; i++) { \
5177 M[i][0].path = M_MATRIX; \
5178 for (j = 1; j <= nB; j++) { \
5179 /* Remove traces to unreachable points. */ \
5180 trace = M[i][j].trace; \
5181 if (!(M[i-1][j-1].path & M_MATRIX)) trace &= ~M_MATRIX; \
5182 if (!(M[i-1][j-1].path & Ix_MATRIX)) trace &= ~Ix_MATRIX; \
5183 if (!(M[i-1][j-1].path & Iy_MATRIX)) trace &= ~Iy_MATRIX; \
5184 if (trace & (STARTPOINT | M_MATRIX | Ix_MATRIX | Iy_MATRIX)) { \
5185 /* The point is reachable. */ \
5186 if (trace & ENDPOINT) M[i][j].path = 0; /* no extensions after ENDPOINT */ \
5187 else M[i][j].path |= M_MATRIX; \
5188 } \
5189 else { \
5190 /* The point is not reachable. Then it is not a STARTPOINT, \
5191 * all traces from it can be removed, and it cannot act as \
5192 * an ENDPOINT. */ \
5193 M[i][j].path &= ~M_MATRIX; \
5194 trace = 0; \
5195 } \
5196 M[i][j].trace = trace; \
5197 trace = gaps[i][j].Ix; \
5198 if (!(M[i-1][j].path & M_MATRIX)) trace &= ~M_MATRIX; \
5199 if (!(M[i-1][j].path & Ix_MATRIX)) trace &= ~Ix_MATRIX; \
5200 if (!(M[i-1][j].path & Iy_MATRIX)) trace &= ~Iy_MATRIX; \
5201 if (trace & (M_MATRIX | Ix_MATRIX | Iy_MATRIX)) { \
5202 /* The point is reachable. */ \
5203 M[i][j].path |= Ix_MATRIX; \
5204 } \
5205 else { \
5206 /* The point is not reachable. Then \
5207 * all traces from it can be removed. */ \
5208 M[i][j].path &= ~Ix_MATRIX; \
5209 trace = 0; \
5210 } \
5211 gaps[i][j].Ix = trace; \
5212 trace = gaps[i][j].Iy; \
5213 if (!(M[i][j-1].path & M_MATRIX)) trace &= ~M_MATRIX; \
5214 if (!(M[i][j-1].path & Ix_MATRIX)) trace &= ~Ix_MATRIX; \
5215 if (!(M[i][j-1].path & Iy_MATRIX)) trace &= ~Iy_MATRIX; \
5216 if (trace & (M_MATRIX | Ix_MATRIX | Iy_MATRIX)) { \
5217 /* The point is reachable. */ \
5218 M[i][j].path |= Iy_MATRIX; \
5219 } \
5220 else { \
5221 /* The point is not reachable. Then \
5222 * all traces from it can be removed. */ \
5223 M[i][j].path &= ~Iy_MATRIX; \
5224 trace = 0; \
5225 } \
5226 gaps[i][j].Iy = trace; \
5227 } \
5228 } \
5229 \
5230 /* traceback */ \
5231 if (maximum == 0) M[0][0].path = DONE; \
5232 else M[0][0].path = 0; \
5233 return Py_BuildValue("fN", maximum, paths); \
5234 \
5235 exit: \
5236 Py_DECREF(paths); \
5237 if (M_row) PyMem_Free(M_row); \
5238 if (Ix_row) PyMem_Free(Ix_row); \
5239 if (Iy_row) PyMem_Free(Iy_row); \
5240 return PyErr_NoMemory(); \
5241
5242
5243 #define WATERMANSMITHBEYER_ENTER_SCORE \
5244 int i; \
5245 int j = 0; \
5246 int k; \
5247 int kA; \
5248 int kB; \
5249 double** M = NULL; \
5250 double** Ix = NULL; \
5251 double** Iy = NULL; \
5252 double score = 0.0; \
5253 double gapscore = 0.0; \
5254 double temp; \
5255 int ok = 1; \
5256 PyObject* result = NULL; \
5257 \
5258 /* Waterman-Smith-Beyer algorithm */ \
5259 M = PyMem_Malloc((nA+1)*sizeof(double*)); \
5260 if (!M) goto exit; \
5261 Ix = PyMem_Malloc((nA+1)*sizeof(double*)); \
5262 if (!Ix) goto exit; \
5263 Iy = PyMem_Malloc((nA+1)*sizeof(double*)); \
5264 if (!Iy) goto exit; \
5265 for (i = 0; i <= nA; i++) { \
5266 M[i] = PyMem_Malloc((nB+1)*sizeof(double)); \
5267 if (!M[i]) goto exit; \
5268 Ix[i] = PyMem_Malloc((nB+1)*sizeof(double)); \
5269 if (!Ix[i]) goto exit; \
5270 Iy[i] = PyMem_Malloc((nB+1)*sizeof(double)); \
5271 if (!Iy[i]) goto exit; \
5272 } \
5273
5274
5275 #define WATERMANSMITHBEYER_GLOBAL_SCORE(align_score, query_gap_start) \
5276 /* The top row of the score matrix is a special case, \
5277 * as there are no previously aligned characters. \
5278 */ \
5279 M[0][0] = 0; \
5280 Ix[0][0] = -DBL_MAX; \
5281 Iy[0][0] = -DBL_MAX; \
5282 for (i = 1; i <= nA; i++) { \
5283 M[i][0] = -DBL_MAX; \
5284 Iy[i][0] = -DBL_MAX; \
5285 ok = _call_query_gap_function(self, query_gap_start, i, &score); \
5286 if (!ok) goto exit; \
5287 Ix[i][0] = score; \
5288 } \
5289 for (j = 1; j <= nB; j++) { \
5290 M[0][j] = -DBL_MAX; \
5291 Ix[0][j] = -DBL_MAX; \
5292 ok = _call_target_gap_function(self, 0, j, &score); \
5293 if (!ok) goto exit; \
5294 Iy[0][j] = score; \
5295 } \
5296 for (i = 1; i <= nA; i++) { \
5297 kA = sA[i-1]; \
5298 for (j = 1; j <= nB; j++) { \
5299 kB = sB[j-1]; \
5300 SELECT_SCORE_GLOBAL(M[i-1][j-1], Ix[i-1][j-1], Iy[i-1][j-1]); \
5301 M[i][j] = score + (align_score); \
5302 score = -DBL_MAX; \
5303 for (k = 1; k <= i; k++) { \
5304 ok = _call_query_gap_function(self, query_gap_start, k, &gapscore); \
5305 if (!ok) goto exit; \
5306 SELECT_SCORE_WATERMAN_SMITH_BEYER(M[i-k][j], Iy[i-k][j]); \
5307 } \
5308 Ix[i][j] = score; \
5309 score = -DBL_MAX; \
5310 for (k = 1; k <= j; k++) { \
5311 ok = _call_target_gap_function(self, i, k, &gapscore); \
5312 if (!ok) goto exit; \
5313 SELECT_SCORE_WATERMAN_SMITH_BEYER(M[i][j-k], Ix[i][j-k]); \
5314 } \
5315 Iy[i][j] = score; \
5316 } \
5317 } \
5318 SELECT_SCORE_GLOBAL(M[nA][nB], Ix[nA][nB], Iy[nA][nB]); \
5319 \
5320 result = PyFloat_FromDouble(score); \
5321
5322
5323 #define WATERMANSMITHBEYER_LOCAL_SCORE(align_score, query_gap_start) \
5324 /* The top row of the score matrix is a special case, \
5325 * as there are no previously aligned characters. \
5326 */ \
5327 M[0][0] = 0; \
5328 Ix[0][0] = -DBL_MAX; \
5329 Iy[0][0] = -DBL_MAX; \
5330 for (i = 1; i <= nA; i++) { \
5331 M[i][0] = -DBL_MAX; \
5332 Ix[i][0] = 0; \
5333 Iy[i][0] = -DBL_MAX; \
5334 } \
5335 for (j = 1; j <= nB; j++) { \
5336 M[0][j] = -DBL_MAX; \
5337 Ix[0][j] = -DBL_MAX; \
5338 Iy[0][j] = 0; \
5339 } \
5340 for (i = 1; i <= nA; i++) { \
5341 kA = sA[i-1]; \
5342 for (j = 1; j <= nB; j++) { \
5343 kB = sB[j-1]; \
5344 SELECT_SCORE_GOTOH_LOCAL_ALIGN(M[i-1][j-1], \
5345 Ix[i-1][j-1], \
5346 Iy[i-1][j-1], \
5347 (align_score)); \
5348 M[i][j] = score; \
5349 if (i == nA || j == nB) { \
5350 Ix[i][j] = 0; \
5351 Iy[i][j] = 0; \
5352 continue; \
5353 } \
5354 score = 0.0; \
5355 for (k = 1; k <= i; k++) { \
5356 ok = _call_query_gap_function(self, query_gap_start, k, &gapscore); \
5357 SELECT_SCORE_WATERMAN_SMITH_BEYER(M[i-k][j], Iy[i-k][j]); \
5358 if (!ok) goto exit; \
5359 } \
5360 if (score > maximum) maximum = score; \
5361 Ix[i][j] = score; \
5362 score = 0.0; \
5363 for (k = 1; k <= j; k++) { \
5364 ok = _call_target_gap_function(self, i, k, &gapscore); \
5365 if (!ok) goto exit; \
5366 SELECT_SCORE_WATERMAN_SMITH_BEYER(M[i][j-k], Ix[i][j-k]); \
5367 } \
5368 if (score > maximum) maximum = score; \
5369 Iy[i][j] = score; \
5370 } \
5371 } \
5372 SELECT_SCORE_GLOBAL(M[nA][nB], Ix[nA][nB], Iy[nA][nB]); \
5373 if (score > maximum) maximum = score; \
5374 result = PyFloat_FromDouble(maximum); \
5375
5376
5377 #define WATERMANSMITHBEYER_EXIT_SCORE \
5378 exit: \
5379 if (M) { \
5380 /* If M is NULL, then Ix is also NULL. */ \
5381 if (Ix) { \
5382 /* If Ix is NULL, then Iy is also NULL. */ \
5383 if (Iy) { \
5384 /* If Iy is NULL, then M[i], Ix[i], and Iy[i] are \
5385 * also NULL. */ \
5386 for (i = 0; i <= nA; i++) { \
5387 if (!M[i]) break; \
5388 PyMem_Free(M[i]); \
5389 if (!Ix[i]) break; \
5390 PyMem_Free(Ix[i]); \
5391 if (!Iy[i]) break; \
5392 PyMem_Free(Iy[i]); \
5393 } \
5394 PyMem_Free(Iy); \
5395 } \
5396 PyMem_Free(Ix); \
5397 } \
5398 PyMem_Free(M); \
5399 } \
5400 if (!ok) return NULL; \
5401 if (!result) return PyErr_NoMemory(); \
5402 return result; \
5403
5404
5405 #define WATERMANSMITHBEYER_ENTER_ALIGN(mode) \
5406 int i; \
5407 int j = 0; \
5408 int gap; \
5409 int kA; \
5410 int kB; \
5411 const double epsilon = self->epsilon; \
5412 Trace** M; \
5413 TraceGapsWatermanSmithBeyer** gaps; \
5414 double** M_row; \
5415 double** Ix_row; \
5416 double** Iy_row; \
5417 int ng; \
5418 int nm; \
5419 double score; \
5420 double gapscore; \
5421 double temp; \
5422 int trace; \
5423 int* gapM; \
5424 int* gapXY; \
5425 int ok = 1; \
5426 PathGenerator* paths = NULL; \
5427 \
5428 /* Waterman-Smith-Beyer algorithm */ \
5429 paths = PathGenerator_create_WSB(nA, nB, mode, strand); \
5430 if (!paths) return NULL; \
5431 M = paths->M; \
5432 gaps = paths->gaps.waterman_smith_beyer; \
5433 M_row = PyMem_Malloc((nA+1)*sizeof(double*)); \
5434 if (!M_row) goto exit; \
5435 Ix_row = PyMem_Malloc((nA+1)*sizeof(double*)); \
5436 if (!Ix_row) goto exit; \
5437 Iy_row = PyMem_Malloc((nA+1)*sizeof(double*)); \
5438 if (!Iy_row) goto exit; \
5439 for (i = 0; i <= nA; i++) { \
5440 M_row[i] = PyMem_Malloc((nB+1)*sizeof(double)); \
5441 if (!M_row[i]) goto exit; \
5442 Ix_row[i] = PyMem_Malloc((nB+1)*sizeof(double)); \
5443 if (!Ix_row[i]) goto exit; \
5444 Iy_row[i] = PyMem_Malloc((nB+1)*sizeof(double)); \
5445 if (!Iy_row[i]) goto exit; \
5446 } \
5447
5448
5449 #define WATERMANSMITHBEYER_GLOBAL_ALIGN(align_score, query_gap_start) \
5450 M_row[0][0] = 0; \
5451 Ix_row[0][0] = -DBL_MAX; \
5452 Iy_row[0][0] = -DBL_MAX; \
5453 for (i = 1; i <= nA; i++) { \
5454 M_row[i][0] = -DBL_MAX; \
5455 Iy_row[i][0] = -DBL_MAX; \
5456 ok = _call_query_gap_function(self, query_gap_start, i, &score); \
5457 if (!ok) goto exit; \
5458 Ix_row[i][0] = score; \
5459 } \
5460 for (j = 1; j <= nB; j++) { \
5461 M_row[0][j] = -DBL_MAX; \
5462 Ix_row[0][j] = -DBL_MAX; \
5463 ok = _call_target_gap_function(self, query_gap_start, j, &score); \
5464 if (!ok) goto exit; \
5465 Iy_row[0][j] = score; \
5466 } \
5467 for (i = 1; i <= nA; i++) { \
5468 kA = sA[i-1]; \
5469 for (j = 1; j <= nB; j++) { \
5470 kB = sB[j-1]; \
5471 SELECT_TRACE_WATERMAN_SMITH_BEYER_GLOBAL_ALIGN((align_score)); \
5472 gapM = PyMem_Malloc((i+1)*sizeof(int)); \
5473 if (!gapM) goto exit; \
5474 gaps[i][j].MIx = gapM; \
5475 gapXY = PyMem_Malloc((i+1)*sizeof(int)); \
5476 if (!gapXY) goto exit; \
5477 gaps[i][j].IyIx = gapXY; \
5478 nm = 0; \
5479 ng = 0; \
5480 score = -DBL_MAX; \
5481 for (gap = 1; gap <= i; gap++) { \
5482 ok = _call_query_gap_function(self, query_gap_start, gap, &gapscore); \
5483 if (!ok) goto exit; \
5484 SELECT_TRACE_WATERMAN_SMITH_BEYER_GAP(M_row[i-gap][j], \
5485 Iy_row[i-gap][j]); \
5486 } \
5487 gapM = PyMem_Realloc(gapM, (nm+1)*sizeof(int)); \
5488 if (!gapM) goto exit; \
5489 gaps[i][j].MIx = gapM; \
5490 gapM[nm] = 0; \
5491 gapXY = PyMem_Realloc(gapXY, (ng+1)*sizeof(int)); \
5492 if (!gapXY) goto exit; \
5493 gapXY[ng] = 0; \
5494 gaps[i][j].IyIx = gapXY; \
5495 Ix_row[i][j] = score; \
5496 gapM = PyMem_Malloc((j+1)*sizeof(int)); \
5497 if (!gapM) goto exit; \
5498 gaps[i][j].MIy = gapM; \
5499 gapXY = PyMem_Malloc((j+1)*sizeof(int)); \
5500 if (!gapXY) goto exit; \
5501 gaps[i][j].IxIy = gapXY; \
5502 nm = 0; \
5503 ng = 0; \
5504 score = -DBL_MAX; \
5505 for (gap = 1; gap <= j; gap++) { \
5506 ok = _call_target_gap_function(self, i, gap, &gapscore); \
5507 if (!ok) goto exit; \
5508 SELECT_TRACE_WATERMAN_SMITH_BEYER_GAP(M_row[i][j-gap], \
5509 Ix_row[i][j-gap]); \
5510 } \
5511 Iy_row[i][j] = score; \
5512 gapM = PyMem_Realloc(gapM, (nm+1)*sizeof(int)); \
5513 if (!gapM) goto exit; \
5514 gaps[i][j].MIy = gapM; \
5515 gapM[nm] = 0; \
5516 gapXY = PyMem_Realloc(gapXY, (ng+1)*sizeof(int)); \
5517 if (!gapXY) goto exit; \
5518 gaps[i][j].IxIy = gapXY; \
5519 gapXY[ng] = 0; \
5520 } \
5521 } \
5522 /* traceback */ \
5523 SELECT_SCORE_GLOBAL(M_row[nA][nB], Ix_row[nA][nB], Iy_row[nA][nB]); \
5524 M[nA][nB].path = 0; \
5525 if (M_row[nA][nB] < score - epsilon) M[nA][nB].trace = 0; \
5526 if (Ix_row[nA][nB] < score - epsilon) { \
5527 gapM = PyMem_Realloc(gaps[nA][nB].MIx, sizeof(int)); \
5528 if (!gapM) goto exit; \
5529 gapM[0] = 0; \
5530 gaps[nA][nB].MIx = gapM; \
5531 gapXY = PyMem_Realloc(gaps[nA][nB].IyIx, sizeof(int)); \
5532 if (!gapXY) goto exit; \
5533 gapXY[0] = 0; \
5534 gaps[nA][nB].IyIx = gapXY; \
5535 } \
5536 if (Iy_row[nA][nB] < score - epsilon) { \
5537 gapM = PyMem_Realloc(gaps[nA][nB].MIy, sizeof(int)); \
5538 if (!gapM) goto exit; \
5539 gapM[0] = 0; \
5540 gaps[nA][nB].MIy = gapM; \
5541 gapXY = PyMem_Realloc(gaps[nA][nB].IxIy, sizeof(int)); \
5542 if (!gapXY) goto exit; \
5543 gapXY[0] = 0; \
5544 gaps[nA][nB].IxIy = gapXY; \
5545 } \
5546 for (i = 0; i <= nA; i++) { \
5547 PyMem_Free(M_row[i]); \
5548 PyMem_Free(Ix_row[i]); \
5549 PyMem_Free(Iy_row[i]); \
5550 } \
5551 PyMem_Free(M_row); \
5552 PyMem_Free(Ix_row); \
5553 PyMem_Free(Iy_row); \
5554 return Py_BuildValue("fN", score, paths); \
5555
5556
5557 #define WATERMANSMITHBEYER_LOCAL_ALIGN(align_score, query_gap_start) \
5558 M_row[0][0] = 0; \
5559 Ix_row[0][0] = -DBL_MAX; \
5560 Iy_row[0][0] = -DBL_MAX; \
5561 for (i = 1; i <= nA; i++) { \
5562 M_row[i][0] = 0; \
5563 Ix_row[i][0] = -DBL_MAX; \
5564 Iy_row[i][0] = -DBL_MAX; \
5565 } \
5566 for (i = 1; i <= nB; i++) { \
5567 M_row[0][i] = 0; \
5568 Ix_row[0][i] = -DBL_MAX; \
5569 Iy_row[0][i] = -DBL_MAX; \
5570 } \
5571 for (i = 1; i <= nA; i++) { \
5572 kA = sA[i-1]; \
5573 for (j = 1; j <= nB; j++) { \
5574 kB = sB[j-1]; \
5575 nm = 0; \
5576 ng = 0; \
5577 SELECT_TRACE_WATERMAN_SMITH_BEYER_ALIGN( \
5578 M_row[i-1][j-1], \
5579 Ix_row[i-1][j-1], \
5580 Iy_row[i-1][j-1], \
5581 (align_score)); \
5582 M[i][j].path = 0; \
5583 if (i == nA || j == nB) { \
5584 Ix_row[i][j] = score; \
5585 gaps[i][j].MIx = NULL; \
5586 gaps[i][j].IyIx = NULL; \
5587 gaps[i][j].MIy = NULL; \
5588 gaps[i][j].IxIy = NULL; \
5589 Iy_row[i][j] = score; \
5590 continue; \
5591 } \
5592 gapM = PyMem_Malloc((i+1)*sizeof(int)); \
5593 if (!gapM) goto exit; \
5594 gaps[i][j].MIx = gapM; \
5595 gapXY = PyMem_Malloc((i+1)*sizeof(int)); \
5596 if (!gapXY) goto exit; \
5597 gaps[i][j].IyIx = gapXY; \
5598 score = -DBL_MAX; \
5599 for (gap = 1; gap <= i; gap++) { \
5600 ok = _call_query_gap_function(self, query_gap_start, gap, &gapscore); \
5601 if (!ok) goto exit; \
5602 SELECT_TRACE_WATERMAN_SMITH_BEYER_GAP(M_row[i-gap][j], \
5603 Iy_row[i-gap][j]); \
5604 } \
5605 if (score < epsilon) { \
5606 score = -DBL_MAX; \
5607 nm = 0; \
5608 ng = 0; \
5609 } \
5610 else if (score > maximum) maximum = score; \
5611 gapM[nm] = 0; \
5612 gapXY[ng] = 0; \
5613 Ix_row[i][j] = score; \
5614 M[i][j].path = 0; \
5615 gapM = PyMem_Realloc(gapM, (nm+1)*sizeof(int)); \
5616 if (!gapM) goto exit; \
5617 gaps[i][j].MIx = gapM; \
5618 gapM[nm] = 0; \
5619 gapXY = PyMem_Realloc(gapXY, (ng+1)*sizeof(int)); \
5620 if (!gapXY) goto exit; \
5621 gaps[i][j].IyIx = gapXY; \
5622 gapXY[ng] = 0; \
5623 gapM = PyMem_Malloc((j+1)*sizeof(int)); \
5624 if (!gapM) goto exit; \
5625 gaps[i][j].MIy = gapM; \
5626 gapXY = PyMem_Malloc((j+1)*sizeof(int)); \
5627 if (!gapXY) goto exit; \
5628 gaps[i][j].IxIy = gapXY; \
5629 nm = 0; \
5630 ng = 0; \
5631 score = -DBL_MAX; \
5632 gapM[0] = 0; \
5633 for (gap = 1; gap <= j; gap++) { \
5634 ok = _call_target_gap_function(self, i, gap, &gapscore); \
5635 if (!ok) goto exit; \
5636 SELECT_TRACE_WATERMAN_SMITH_BEYER_GAP(M_row[i][j-gap], \
5637 Ix_row[i][j-gap]); \
5638 } \
5639 if (score < epsilon) { \
5640 score = -DBL_MAX; \
5641 nm = 0; \
5642 ng = 0; \
5643 } \
5644 else if (score > maximum) maximum = score; \
5645 gapM = PyMem_Realloc(gapM, (nm+1)*sizeof(int)); \
5646 if (!gapM) goto exit; \
5647 gaps[i][j].MIy = gapM; \
5648 gapXY = PyMem_Realloc(gapXY, (ng+1)*sizeof(int)); \
5649 if (!gapXY) goto exit; \
5650 gaps[i][j].IxIy = gapXY; \
5651 gapM[nm] = 0; \
5652 gapXY[ng] = 0; \
5653 Iy_row[i][j] = score; \
5654 M[i][j].path = 0; \
5655 } \
5656 } \
5657 for (i = 0; i <= nA; i++) PyMem_Free(M_row[i]); \
5658 PyMem_Free(M_row); \
5659 for (i = 0; i <= nA; i++) PyMem_Free(Ix_row[i]); \
5660 PyMem_Free(Ix_row); \
5661 for (i = 0; i <= nA; i++) PyMem_Free(Iy_row[i]); \
5662 PyMem_Free(Iy_row); \
5663 \
5664 /* As we don't allow zero-score extensions to alignments, \
5665 * we need to remove all traces towards an ENDPOINT. \
5666 * In addition, some points then won't have any path to a STARTPOINT. \
5667 * Here, use path as a temporary variable to indicate if the point \
5668 * is reachable from a STARTPOINT. If it is unreachable, remove all \
5669 * traces from it, and don't allow it to be an ENDPOINT. It may still \
5670 * be a valid STARTPOINT. */ \
5671 for (j = 0; j <= nB; j++) M[0][j].path = M_MATRIX; \
5672 for (i = 1; i <= nA; i++) { \
5673 M[i][0].path = M_MATRIX; \
5674 for (j = 1; j <= nB; j++) { \
5675 /* Remove traces to unreachable points. */ \
5676 trace = M[i][j].trace; \
5677 if (!(M[i-1][j-1].path & M_MATRIX)) trace &= ~M_MATRIX; \
5678 if (!(M[i-1][j-1].path & Ix_MATRIX)) trace &= ~Ix_MATRIX; \
5679 if (!(M[i-1][j-1].path & Iy_MATRIX)) trace &= ~Iy_MATRIX; \
5680 if (trace & (STARTPOINT | M_MATRIX | Ix_MATRIX | Iy_MATRIX)) { \
5681 /* The point is reachable. */ \
5682 if (trace & ENDPOINT) M[i][j].path = 0; /* no extensions after ENDPOINT */ \
5683 else M[i][j].path |= M_MATRIX; \
5684 } \
5685 else { \
5686 /* The point is not reachable. Then it is not a STARTPOINT, \
5687 * all traces from it can be removed, and it cannot act as \
5688 * an ENDPOINT. */ \
5689 M[i][j].path &= ~M_MATRIX; \
5690 trace = 0; \
5691 } \
5692 M[i][j].trace = trace; \
5693 if (i == nA || j == nB) continue; \
5694 gapM = gaps[i][j].MIx; \
5695 gapXY = gaps[i][j].IyIx; \
5696 nm = 0; \
5697 ng = 0; \
5698 for (im = 0; (gap = gapM[im]); im++) \
5699 if (M[i-gap][j].path & M_MATRIX) gapM[nm++] = gap; \
5700 gapM = PyMem_Realloc(gapM, (nm+1)*sizeof(int)); \
5701 if (!gapM) goto exit; \
5702 gapM[nm] = 0; \
5703 gaps[i][j].MIx = gapM; \
5704 for (im = 0; (gap = gapXY[im]); im++) \
5705 if (M[i-gap][j].path & Iy_MATRIX) gapXY[ng++] = gap; \
5706 gapXY = PyMem_Realloc(gapXY, (ng+1)*sizeof(int)); \
5707 if (!gapXY) goto exit; \
5708 gapXY[ng] = 0; \
5709 gaps[i][j].IyIx = gapXY; \
5710 if (nm==0 && ng==0) M[i][j].path &= ~Ix_MATRIX; /* not reachable */ \
5711 else M[i][j].path |= Ix_MATRIX; /* reachable */ \
5712 gapM = gaps[i][j].MIy; \
5713 gapXY = gaps[i][j].IxIy; \
5714 nm = 0; \
5715 ng = 0; \
5716 for (im = 0; (gap = gapM[im]); im++) \
5717 if (M[i][j-gap].path & M_MATRIX) gapM[nm++] = gap; \
5718 gapM = PyMem_Realloc(gapM, (nm+1)*sizeof(int)); \
5719 if (!gapM) goto exit; \
5720 gapM[nm] = 0; \
5721 gaps[i][j].MIy = gapM; \
5722 for (im = 0; (gap = gapXY[im]); im++) \
5723 if (M[i][j-gap].path & Ix_MATRIX) gapXY[ng++] = gap; \
5724 gapXY = PyMem_Realloc(gapXY, (ng+1)*sizeof(int)); \
5725 if (!gapXY) goto exit; \
5726 gapXY[ng] = 0; \
5727 gaps[i][j].IxIy = gapXY; \
5728 if (nm==0 && ng==0) M[i][j].path &= ~Iy_MATRIX; /* not reachable */ \
5729 else M[i][j].path |= Iy_MATRIX; /* reachable */ \
5730 } \
5731 } \
5732 /* traceback */ \
5733 if (maximum == 0) M[0][0].path = DONE; \
5734 else M[0][0].path = 0; \
5735 return Py_BuildValue("fN", maximum, paths); \
5736
5737
5738 #define WATERMANSMITHBEYER_EXIT_ALIGN \
5739 exit: \
5740 if (ok) /* otherwise, an exception was already set */ \
5741 PyErr_SetNone(PyExc_MemoryError); \
5742 Py_DECREF(paths); \
5743 if (M_row) { \
5744 /* If M is NULL, then Ix is also NULL. */ \
5745 if (Ix_row) { \
5746 /* If Ix is NULL, then Iy is also NULL. */ \
5747 if (Iy_row) { \
5748 /* If Iy is NULL, then M[i], Ix[i], and Iy[i] are also NULL. */ \
5749 for (i = 0; i <= nA; i++) { \
5750 if (!M_row[i]) break; \
5751 PyMem_Free(M_row[i]); \
5752 if (!Ix_row[i]) break; \
5753 PyMem_Free(Ix_row[i]); \
5754 if (!Iy_row[i]) break; \
5755 PyMem_Free(Iy_row[i]); \
5756 } \
5757 PyMem_Free(Iy_row); \
5758 } \
5759 PyMem_Free(Ix_row); \
5760 } \
5761 PyMem_Free(M_row); \
5762 } \
5763 return NULL; \
5764
5765
5766 /* -------------- allocation & deallocation ------------- */
5767
5768 static PathGenerator*
PathGenerator_create_NWSW(Py_ssize_t nA,Py_ssize_t nB,Mode mode,unsigned char strand)5769 PathGenerator_create_NWSW(Py_ssize_t nA, Py_ssize_t nB, Mode mode, unsigned char strand)
5770 {
5771 int i;
5772 unsigned char trace = 0;
5773 Trace** M;
5774 PathGenerator* paths;
5775
5776 paths = (PathGenerator*)PyType_GenericAlloc(&PathGenerator_Type, 0);
5777 if (!paths) return NULL;
5778
5779 paths->iA = 0;
5780 paths->iB = 0;
5781 paths->nA = nA;
5782 paths->nB = nB;
5783 paths->M = NULL;
5784 paths->gaps.gotoh = NULL;
5785 paths->gaps.waterman_smith_beyer = NULL;
5786 paths->algorithm = NeedlemanWunschSmithWaterman;
5787 paths->mode = mode;
5788 paths->length = 0;
5789 paths->strand = strand;
5790
5791 M = PyMem_Malloc((nA+1)*sizeof(Trace*));
5792 paths->M = M;
5793 if (!M) goto exit;
5794 switch (mode) {
5795 case Global: trace = VERTICAL; break;
5796 case Local: trace = STARTPOINT; break;
5797 }
5798 for (i = 0; i <= nA; i++) {
5799 M[i] = PyMem_Malloc((nB+1)*sizeof(Trace));
5800 if (!M[i]) goto exit;
5801 M[i][0].trace = trace;
5802 }
5803 if (mode == Global) {
5804 M[0][0].trace = 0;
5805 trace = HORIZONTAL;
5806 }
5807 for (i = 1; i <= nB; i++) M[0][i].trace = trace;
5808 M[0][0].path = 0;
5809 return paths;
5810 exit:
5811 Py_DECREF(paths);
5812 PyErr_SetNone(PyExc_MemoryError);
5813 return NULL;
5814 }
5815
5816 static PathGenerator*
PathGenerator_create_Gotoh(Py_ssize_t nA,Py_ssize_t nB,Mode mode,unsigned char strand)5817 PathGenerator_create_Gotoh(Py_ssize_t nA, Py_ssize_t nB, Mode mode, unsigned char strand)
5818 {
5819 int i;
5820 unsigned char trace;
5821 Trace** M;
5822 TraceGapsGotoh** gaps;
5823 PathGenerator* paths;
5824
5825 switch (mode) {
5826 case Global: trace = 0; break;
5827 case Local: trace = STARTPOINT; break;
5828 default:
5829 /* Should not happen, but the compiler has no way of knowing that,
5830 * as the enum Mode does not restrict the value of mode, which can
5831 * be any integer. Include default: here to prevent compiler
5832 * warnings.
5833 */
5834 PyErr_Format(PyExc_RuntimeError,
5835 "mode has unexpected value %d", mode);
5836 return NULL;
5837 }
5838
5839 paths = (PathGenerator*)PyType_GenericAlloc(&PathGenerator_Type, 0);
5840 if (!paths) return NULL;
5841
5842 paths->iA = 0;
5843 paths->iB = 0;
5844 paths->nA = nA;
5845 paths->nB = nB;
5846 paths->M = NULL;
5847 paths->gaps.gotoh = NULL;
5848 paths->algorithm = Gotoh;
5849 paths->mode = mode;
5850 paths->length = 0;
5851 paths->strand = strand;
5852
5853 M = PyMem_Malloc((nA+1)*sizeof(Trace*));
5854 if (!M) goto exit;
5855 paths->M = M;
5856 for (i = 0; i <= nA; i++) {
5857 M[i] = PyMem_Malloc((nB+1)*sizeof(Trace));
5858 if (!M[i]) goto exit;
5859 M[i][0].trace = trace;
5860 }
5861 gaps = PyMem_Malloc((nA+1)*sizeof(TraceGapsGotoh*));
5862 if (!gaps) goto exit;
5863 paths->gaps.gotoh = gaps;
5864 for (i = 0; i <= nA; i++) {
5865 gaps[i] = PyMem_Malloc((nB+1)*sizeof(TraceGapsGotoh));
5866 if (!gaps[i]) goto exit;
5867 }
5868
5869 gaps[0][0].Ix = 0;
5870 gaps[0][0].Iy = 0;
5871 if (mode == Global) {
5872 for (i = 1; i <= nA; i++) {
5873 gaps[i][0].Ix = Ix_MATRIX;
5874 gaps[i][0].Iy = 0;
5875 }
5876 gaps[1][0].Ix = M_MATRIX;
5877 for (i = 1; i <= nB; i++) {
5878 M[0][i].trace = 0;
5879 gaps[0][i].Ix = 0;
5880 gaps[0][i].Iy = Iy_MATRIX;
5881 }
5882 gaps[0][1].Iy = M_MATRIX;
5883 }
5884 else if (mode == Local) {
5885 for (i = 1; i < nA; i++) {
5886 gaps[i][0].Ix = 0;
5887 gaps[i][0].Iy = 0;
5888 }
5889 for (i = 1; i <= nB; i++) {
5890 M[0][i].trace = trace;
5891 gaps[0][i].Ix = 0;
5892 gaps[0][i].Iy = 0;
5893 }
5894 }
5895 M[0][0].path = 0;
5896
5897 return paths;
5898 exit:
5899 Py_DECREF(paths);
5900 PyErr_SetNone(PyExc_MemoryError);
5901 return NULL;
5902 }
5903
5904 static PathGenerator*
PathGenerator_create_WSB(Py_ssize_t nA,Py_ssize_t nB,Mode mode,unsigned char strand)5905 PathGenerator_create_WSB(Py_ssize_t nA, Py_ssize_t nB, Mode mode, unsigned char strand)
5906 {
5907 int i, j;
5908 int* trace;
5909 Trace** M = NULL;
5910 TraceGapsWatermanSmithBeyer** gaps = NULL;
5911 PathGenerator* paths;
5912
5913 paths = (PathGenerator*)PyType_GenericAlloc(&PathGenerator_Type, 0);
5914 if (!paths) return NULL;
5915
5916 paths->iA = 0;
5917 paths->iB = 0;
5918 paths->nA = nA;
5919 paths->nB = nB;
5920 paths->M = NULL;
5921 paths->gaps.waterman_smith_beyer = NULL;
5922 paths->algorithm = WatermanSmithBeyer;
5923 paths->mode = mode;
5924 paths->length = 0;
5925 paths->strand = strand;
5926
5927 M = PyMem_Malloc((nA+1)*sizeof(Trace*));
5928 if (!M) goto exit;
5929 paths->M = M;
5930 for (i = 0; i <= nA; i++) {
5931 M[i] = PyMem_Malloc((nB+1)*sizeof(Trace));
5932 if (!M[i]) goto exit;
5933 }
5934 gaps = PyMem_Malloc((nA+1)*sizeof(TraceGapsWatermanSmithBeyer*));
5935 if (!gaps) goto exit;
5936 paths->gaps.waterman_smith_beyer = gaps;
5937 for (i = 0; i <= nA; i++) gaps[i] = NULL;
5938 for (i = 0; i <= nA; i++) {
5939 gaps[i] = PyMem_Malloc((nB+1)*sizeof(TraceGapsWatermanSmithBeyer));
5940 if (!gaps[i]) goto exit;
5941 for (j = 0; j <= nB; j++) {
5942 gaps[i][j].MIx = NULL;
5943 gaps[i][j].IyIx = NULL;
5944 gaps[i][j].MIy = NULL;
5945 gaps[i][j].IxIy = NULL;
5946 }
5947 M[i][0].path = 0;
5948 switch (mode) {
5949 case Global:
5950 M[i][0].trace = 0;
5951 trace = PyMem_Malloc(2*sizeof(int));
5952 if (!trace) goto exit;
5953 gaps[i][0].MIx = trace;
5954 trace[0] = i;
5955 trace[1] = 0;
5956 trace = PyMem_Malloc(sizeof(int));
5957 if (!trace) goto exit;
5958 gaps[i][0].IyIx = trace;
5959 trace[0] = 0;
5960 break;
5961 case Local:
5962 M[i][0].trace = STARTPOINT;
5963 break;
5964 }
5965 }
5966 for (i = 1; i <= nB; i++) {
5967 switch (mode) {
5968 case Global:
5969 M[0][i].trace = 0;
5970 trace = PyMem_Malloc(2*sizeof(int));
5971 if (!trace) goto exit;
5972 gaps[0][i].MIy = trace;
5973 trace[0] = i;
5974 trace[1] = 0;
5975 trace = PyMem_Malloc(sizeof(int));
5976 if (!trace) goto exit;
5977 gaps[0][i].IxIy = trace;
5978 trace[0] = 0;
5979 break;
5980 case Local:
5981 M[0][i].trace = STARTPOINT;
5982 break;
5983 }
5984 }
5985 M[0][0].path = 0;
5986 return paths;
5987 exit:
5988 Py_DECREF(paths);
5989 PyErr_SetNone(PyExc_MemoryError);
5990 return NULL;
5991 }
5992
5993 /* ----------------- alignment algorithms ----------------- */
5994
5995 #define MATRIX_SCORE scores[kA*n+kB]
5996 #define COMPARE_SCORE (kA == wildcard || kB == wildcard) ? 0 : (kA == kB) ? match : mismatch
5997
5998
5999 static PyObject*
Aligner_needlemanwunsch_score_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6000 Aligner_needlemanwunsch_score_compare(Aligner* self,
6001 const int* sA, Py_ssize_t nA,
6002 const int* sB, Py_ssize_t nB,
6003 unsigned char strand)
6004 {
6005 const double match = self->match;
6006 const double mismatch = self->mismatch;
6007 const int wildcard = self->wildcard;
6008 NEEDLEMANWUNSCH_SCORE(COMPARE_SCORE);
6009 }
6010
6011 static PyObject*
Aligner_needlemanwunsch_score_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6012 Aligner_needlemanwunsch_score_matrix(Aligner* self,
6013 const int* sA, Py_ssize_t nA,
6014 const int* sB, Py_ssize_t nB,
6015 unsigned char strand)
6016 {
6017 const Py_ssize_t n = self->substitution_matrix.shape[0];
6018 const double* scores = self->substitution_matrix.buf;
6019 NEEDLEMANWUNSCH_SCORE(MATRIX_SCORE);
6020 }
6021
6022 static PyObject*
Aligner_smithwaterman_score_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB)6023 Aligner_smithwaterman_score_compare(Aligner* self,
6024 const int* sA, Py_ssize_t nA,
6025 const int* sB, Py_ssize_t nB)
6026 {
6027 const double match = self->match;
6028 const double mismatch = self->mismatch;
6029 const int wildcard = self->wildcard;
6030 SMITHWATERMAN_SCORE(COMPARE_SCORE);
6031 }
6032
6033 static PyObject*
Aligner_smithwaterman_score_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB)6034 Aligner_smithwaterman_score_matrix(Aligner* self,
6035 const int* sA, Py_ssize_t nA,
6036 const int* sB, Py_ssize_t nB)
6037 {
6038 const Py_ssize_t n = self->substitution_matrix.shape[0];
6039 const double* scores = self->substitution_matrix.buf;
6040 SMITHWATERMAN_SCORE(MATRIX_SCORE);
6041 }
6042
6043 static PyObject*
Aligner_needlemanwunsch_align_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6044 Aligner_needlemanwunsch_align_compare(Aligner* self,
6045 const int* sA, Py_ssize_t nA,
6046 const int* sB, Py_ssize_t nB,
6047 unsigned char strand)
6048 {
6049 const double match = self->match;
6050 const double mismatch = self->mismatch;
6051 const int wildcard = self->wildcard;
6052 NEEDLEMANWUNSCH_ALIGN(COMPARE_SCORE);
6053 }
6054
6055 static PyObject*
Aligner_needlemanwunsch_align_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6056 Aligner_needlemanwunsch_align_matrix(Aligner* self,
6057 const int* sA, Py_ssize_t nA,
6058 const int* sB, Py_ssize_t nB,
6059 unsigned char strand)
6060 {
6061 const Py_ssize_t n = self->substitution_matrix.shape[0];
6062 const double* scores = self->substitution_matrix.buf;
6063 NEEDLEMANWUNSCH_ALIGN(MATRIX_SCORE);
6064 }
6065
6066 static PyObject*
Aligner_smithwaterman_align_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6067 Aligner_smithwaterman_align_compare(Aligner* self,
6068 const int* sA, Py_ssize_t nA,
6069 const int* sB, Py_ssize_t nB,
6070 unsigned char strand)
6071 {
6072 const double match = self->match;
6073 const double mismatch = self->mismatch;
6074 const int wildcard = self->wildcard;
6075 SMITHWATERMAN_ALIGN(COMPARE_SCORE);
6076 }
6077
6078 static PyObject*
Aligner_smithwaterman_align_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6079 Aligner_smithwaterman_align_matrix(Aligner* self,
6080 const int* sA, Py_ssize_t nA,
6081 const int* sB, Py_ssize_t nB,
6082 unsigned char strand)
6083 {
6084 const Py_ssize_t n = self->substitution_matrix.shape[0];
6085 const double* scores = self->substitution_matrix.buf;
6086 SMITHWATERMAN_ALIGN(MATRIX_SCORE);
6087 }
6088
6089 static PyObject*
Aligner_gotoh_global_score_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6090 Aligner_gotoh_global_score_compare(Aligner* self,
6091 const int* sA, Py_ssize_t nA,
6092 const int* sB, Py_ssize_t nB,
6093 unsigned char strand)
6094 {
6095 const double match = self->match;
6096 const double mismatch = self->mismatch;
6097 const int wildcard = self->wildcard;
6098 GOTOH_GLOBAL_SCORE(COMPARE_SCORE);
6099 }
6100
6101 static PyObject*
Aligner_gotoh_global_score_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6102 Aligner_gotoh_global_score_matrix(Aligner* self,
6103 const int* sA, Py_ssize_t nA,
6104 const int* sB, Py_ssize_t nB,
6105 unsigned char strand)
6106 {
6107 const Py_ssize_t n = self->substitution_matrix.shape[0];
6108 const double* scores = self->substitution_matrix.buf;
6109 GOTOH_GLOBAL_SCORE(MATRIX_SCORE);
6110 }
6111
6112 static PyObject*
Aligner_gotoh_local_score_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB)6113 Aligner_gotoh_local_score_compare(Aligner* self,
6114 const int* sA, Py_ssize_t nA,
6115 const int* sB, Py_ssize_t nB)
6116 {
6117 const double match = self->match;
6118 const double mismatch = self->mismatch;
6119 const int wildcard = self->wildcard;
6120 GOTOH_LOCAL_SCORE(COMPARE_SCORE);
6121 }
6122
6123 static PyObject*
Aligner_gotoh_local_score_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB)6124 Aligner_gotoh_local_score_matrix(Aligner* self,
6125 const int* sA, Py_ssize_t nA,
6126 const int* sB, Py_ssize_t nB)
6127 {
6128 const Py_ssize_t n = self->substitution_matrix.shape[0];
6129 const double* scores = self->substitution_matrix.buf;
6130 GOTOH_LOCAL_SCORE(MATRIX_SCORE);
6131 }
6132
6133 static PyObject*
Aligner_gotoh_global_align_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6134 Aligner_gotoh_global_align_compare(Aligner* self,
6135 const int* sA, Py_ssize_t nA,
6136 const int* sB, Py_ssize_t nB,
6137 unsigned char strand)
6138 {
6139 const double match = self->match;
6140 const double mismatch = self->mismatch;
6141 const int wildcard = self->wildcard;
6142 GOTOH_GLOBAL_ALIGN(COMPARE_SCORE);
6143 }
6144
6145 static PyObject*
Aligner_gotoh_global_align_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6146 Aligner_gotoh_global_align_matrix(Aligner* self,
6147 const int* sA, Py_ssize_t nA,
6148 const int* sB, Py_ssize_t nB,
6149 unsigned char strand)
6150 {
6151 const Py_ssize_t n = self->substitution_matrix.shape[0];
6152 const double* scores = self->substitution_matrix.buf;
6153 GOTOH_GLOBAL_ALIGN(MATRIX_SCORE);
6154 }
6155
6156 static PyObject*
Aligner_gotoh_local_align_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6157 Aligner_gotoh_local_align_compare(Aligner* self,
6158 const int* sA, Py_ssize_t nA,
6159 const int* sB, Py_ssize_t nB,
6160 unsigned char strand)
6161 {
6162 const double match = self->match;
6163 const double mismatch = self->mismatch;
6164 const int wildcard = self->wildcard;
6165 GOTOH_LOCAL_ALIGN(COMPARE_SCORE);
6166 }
6167
6168 static PyObject*
Aligner_gotoh_local_align_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6169 Aligner_gotoh_local_align_matrix(Aligner* self,
6170 const int* sA, Py_ssize_t nA,
6171 const int* sB, Py_ssize_t nB,
6172 unsigned char strand)
6173 {
6174 const Py_ssize_t n = self->substitution_matrix.shape[0];
6175 const double* scores = self->substitution_matrix.buf;
6176 GOTOH_LOCAL_ALIGN(MATRIX_SCORE);
6177 }
6178
6179 static int
_call_query_gap_function(Aligner * aligner,int i,int j,double * score)6180 _call_query_gap_function(Aligner* aligner, int i, int j, double* score)
6181 {
6182 double value;
6183 PyObject* result;
6184 PyObject* function = aligner->query_gap_function;
6185 if (!function)
6186 value = aligner->query_internal_open_gap_score
6187 + (j-1) * aligner->query_internal_extend_gap_score;
6188 else {
6189 result = PyObject_CallFunction(function, "ii", i, j);
6190 if (result == NULL) return 0;
6191 value = PyFloat_AsDouble(result);
6192 Py_DECREF(result);
6193 if (value == -1.0 && PyErr_Occurred()) return 0;
6194 }
6195 *score = value;
6196 return 1;
6197 }
6198
6199 static int
_call_target_gap_function(Aligner * aligner,int i,int j,double * score)6200 _call_target_gap_function(Aligner* aligner, int i, int j, double* score)
6201 {
6202 double value;
6203 PyObject* result;
6204 PyObject* function = aligner->target_gap_function;
6205 if (!function)
6206 value = aligner->target_internal_open_gap_score
6207 + (j-1) * aligner->target_internal_extend_gap_score;
6208 else {
6209 result = PyObject_CallFunction(function, "ii", i, j);
6210 if (result == NULL) return 0;
6211 value = PyFloat_AsDouble(result);
6212 Py_DECREF(result);
6213 if (value == -1.0 && PyErr_Occurred()) return 0;
6214 }
6215 *score = value;
6216 return 1;
6217 }
6218
6219 static PyObject*
Aligner_watermansmithbeyer_global_score_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6220 Aligner_watermansmithbeyer_global_score_compare(Aligner* self,
6221 const int* sA, Py_ssize_t nA,
6222 const int* sB, Py_ssize_t nB,
6223 unsigned char strand)
6224 {
6225 const double match = self->match;
6226 const double mismatch = self->mismatch;
6227 const int wildcard = self->wildcard;
6228 WATERMANSMITHBEYER_ENTER_SCORE;
6229 switch (strand) {
6230 case '+': {
6231 WATERMANSMITHBEYER_GLOBAL_SCORE(COMPARE_SCORE, j);
6232 break;
6233 }
6234 case '-': {
6235 WATERMANSMITHBEYER_GLOBAL_SCORE(COMPARE_SCORE, nB-j);
6236 break;
6237 }
6238 }
6239 WATERMANSMITHBEYER_EXIT_SCORE;
6240 }
6241
6242 static PyObject*
Aligner_watermansmithbeyer_global_score_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6243 Aligner_watermansmithbeyer_global_score_matrix(Aligner* self,
6244 const int* sA, Py_ssize_t nA,
6245 const int* sB, Py_ssize_t nB,
6246 unsigned char strand)
6247 {
6248 const Py_ssize_t n = self->substitution_matrix.shape[0];
6249 const double* scores = self->substitution_matrix.buf;
6250 WATERMANSMITHBEYER_ENTER_SCORE;
6251 switch (strand) {
6252 case '+':
6253 WATERMANSMITHBEYER_GLOBAL_SCORE(MATRIX_SCORE, j);
6254 break;
6255 case '-':
6256 WATERMANSMITHBEYER_GLOBAL_SCORE(MATRIX_SCORE, nB-j);
6257 break;
6258 }
6259 WATERMANSMITHBEYER_EXIT_SCORE;
6260 }
6261
6262 static PyObject*
Aligner_watermansmithbeyer_local_score_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6263 Aligner_watermansmithbeyer_local_score_compare(Aligner* self,
6264 const int* sA, Py_ssize_t nA,
6265 const int* sB, Py_ssize_t nB,
6266 unsigned char strand)
6267 {
6268 const double match = self->match;
6269 const double mismatch = self->mismatch;
6270 const int wildcard = self->wildcard;
6271 double maximum = 0.0;
6272 WATERMANSMITHBEYER_ENTER_SCORE;
6273 switch (strand) {
6274 case '+': {
6275 WATERMANSMITHBEYER_LOCAL_SCORE(COMPARE_SCORE, j);
6276 break;
6277 }
6278 case '-': {
6279 WATERMANSMITHBEYER_LOCAL_SCORE(COMPARE_SCORE, nB-j);
6280 break;
6281 }
6282 }
6283 WATERMANSMITHBEYER_EXIT_SCORE;
6284 }
6285
6286 static PyObject*
Aligner_watermansmithbeyer_local_score_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6287 Aligner_watermansmithbeyer_local_score_matrix(Aligner* self,
6288 const int* sA, Py_ssize_t nA,
6289 const int* sB, Py_ssize_t nB,
6290 unsigned char strand)
6291 {
6292 const Py_ssize_t n = self->substitution_matrix.shape[0];
6293 const double* scores = self->substitution_matrix.buf;
6294 double maximum = 0.0;
6295 WATERMANSMITHBEYER_ENTER_SCORE;
6296 switch (strand) {
6297 case '+': {
6298 WATERMANSMITHBEYER_LOCAL_SCORE(MATRIX_SCORE, j);
6299 break;
6300 }
6301 case '-': {
6302 WATERMANSMITHBEYER_LOCAL_SCORE(MATRIX_SCORE, nB-j);
6303 break;
6304 }
6305 }
6306 WATERMANSMITHBEYER_EXIT_SCORE;
6307 }
6308
6309 static PyObject*
Aligner_watermansmithbeyer_global_align_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6310 Aligner_watermansmithbeyer_global_align_compare(Aligner* self,
6311 const int* sA, Py_ssize_t nA,
6312 const int* sB, Py_ssize_t nB,
6313 unsigned char strand)
6314 {
6315 const double match = self->match;
6316 const double mismatch = self->mismatch;
6317 const int wildcard = self->wildcard;
6318 WATERMANSMITHBEYER_ENTER_ALIGN(Global);
6319 switch (strand) {
6320 case '+': {
6321 WATERMANSMITHBEYER_GLOBAL_ALIGN(COMPARE_SCORE, j);
6322 break;
6323 }
6324 case '-': {
6325 WATERMANSMITHBEYER_GLOBAL_ALIGN(COMPARE_SCORE, nB-j);
6326 break;
6327 }
6328 }
6329 WATERMANSMITHBEYER_EXIT_ALIGN;
6330 }
6331
6332 static PyObject*
Aligner_watermansmithbeyer_global_align_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6333 Aligner_watermansmithbeyer_global_align_matrix(Aligner* self,
6334 const int* sA, Py_ssize_t nA,
6335 const int* sB, Py_ssize_t nB,
6336 unsigned char strand)
6337 {
6338 const Py_ssize_t n = self->substitution_matrix.shape[0];
6339 const double* scores = self->substitution_matrix.buf;
6340 WATERMANSMITHBEYER_ENTER_ALIGN(Global);
6341 switch (strand) {
6342 case '+': {
6343 WATERMANSMITHBEYER_GLOBAL_ALIGN(MATRIX_SCORE, j);
6344 break;
6345 }
6346 case '-': {
6347 WATERMANSMITHBEYER_GLOBAL_ALIGN(MATRIX_SCORE, nB-j);
6348 break;
6349 }
6350 }
6351 WATERMANSMITHBEYER_EXIT_ALIGN;
6352 }
6353
6354 static PyObject*
Aligner_watermansmithbeyer_local_align_compare(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6355 Aligner_watermansmithbeyer_local_align_compare(Aligner* self,
6356 const int* sA, Py_ssize_t nA,
6357 const int* sB, Py_ssize_t nB,
6358 unsigned char strand)
6359 {
6360 const double match = self->match;
6361 const double mismatch = self->mismatch;
6362 const int wildcard = self->wildcard;
6363 int im = nA;
6364 int jm = nB;
6365 double maximum = 0;
6366 WATERMANSMITHBEYER_ENTER_ALIGN(Local);
6367 switch (strand) {
6368 case '+': {
6369 WATERMANSMITHBEYER_LOCAL_ALIGN(COMPARE_SCORE, j);
6370 break;
6371 }
6372 case '-': {
6373 WATERMANSMITHBEYER_LOCAL_ALIGN(COMPARE_SCORE, nB-j);
6374 break;
6375 }
6376 }
6377 WATERMANSMITHBEYER_EXIT_ALIGN;
6378 }
6379
6380 static PyObject*
Aligner_watermansmithbeyer_local_align_matrix(Aligner * self,const int * sA,Py_ssize_t nA,const int * sB,Py_ssize_t nB,unsigned char strand)6381 Aligner_watermansmithbeyer_local_align_matrix(Aligner* self,
6382 const int* sA, Py_ssize_t nA,
6383 const int* sB, Py_ssize_t nB,
6384 unsigned char strand)
6385 {
6386 const Py_ssize_t n = self->substitution_matrix.shape[0];
6387 const double* scores = self->substitution_matrix.buf;
6388 int im = nA;
6389 int jm = nB;
6390 double maximum = 0;
6391 WATERMANSMITHBEYER_ENTER_ALIGN(Local);
6392 switch (strand) {
6393 case '+': {
6394 WATERMANSMITHBEYER_LOCAL_ALIGN(MATRIX_SCORE, j);
6395 break;
6396 }
6397 case '-': {
6398 WATERMANSMITHBEYER_LOCAL_ALIGN(MATRIX_SCORE, nB-j);
6399 break;
6400 }
6401 }
6402 WATERMANSMITHBEYER_EXIT_ALIGN;
6403 }
6404
6405 static int*
convert_1bytes_to_ints(const int mapping[],Py_ssize_t n,const unsigned char s[])6406 convert_1bytes_to_ints(const int mapping[], Py_ssize_t n, const unsigned char s[])
6407 {
6408 unsigned char c;
6409 Py_ssize_t i;
6410 int index;
6411 int* indices;
6412 if (n == 0) {
6413 PyErr_SetString(PyExc_ValueError, "sequence has zero length");
6414 return NULL;
6415 }
6416 indices = PyMem_Malloc(n*sizeof(int));
6417 if (!indices) {
6418 PyErr_NoMemory();
6419 return NULL;
6420 }
6421 if (!mapping) for (i = 0; i < n; i++) indices[i] = s[i];
6422 else {
6423 for (i = 0; i < n; i++) {
6424 c = s[i];
6425 index = mapping[(int)c];
6426 if (index == MISSING_LETTER) {
6427 PyErr_SetString(PyExc_ValueError,
6428 "sequence contains letters not in the alphabet");
6429 PyMem_Free(indices);
6430 return NULL;
6431 }
6432 indices[i] = index;
6433 }
6434 }
6435 return indices;
6436 }
6437
6438 static int*
convert_2bytes_to_ints(const int mapping[],Py_ssize_t n,const Py_UCS2 s[])6439 convert_2bytes_to_ints(const int mapping[], Py_ssize_t n, const Py_UCS2 s[])
6440 {
6441 unsigned char c;
6442 Py_ssize_t i;
6443 int index;
6444 int* indices;
6445 if (n == 0) {
6446 PyErr_SetString(PyExc_ValueError, "sequence has zero length");
6447 return NULL;
6448 }
6449 indices = PyMem_Malloc(n*sizeof(int));
6450 if (!indices) {
6451 PyErr_NoMemory();
6452 return NULL;
6453 }
6454 if (!mapping) for (i = 0; i < n; i++) indices[i] = s[i];
6455 else {
6456 for (i = 0; i < n; i++) {
6457 c = s[i];
6458 index = mapping[(int)c];
6459 if (index == MISSING_LETTER) {
6460 PyErr_SetString(PyExc_ValueError,
6461 "sequence contains letters not in the alphabet");
6462 PyMem_Free(indices);
6463 return NULL;
6464 }
6465 indices[i] = index;
6466 }
6467 }
6468 return indices;
6469 }
6470
6471 static int*
convert_4bytes_to_ints(const int mapping[],Py_ssize_t n,const Py_UCS4 s[])6472 convert_4bytes_to_ints(const int mapping[], Py_ssize_t n, const Py_UCS4 s[])
6473 {
6474 unsigned char c;
6475 Py_ssize_t i;
6476 int index;
6477 int* indices;
6478 if (n == 0) {
6479 PyErr_SetString(PyExc_ValueError, "sequence has zero length");
6480 return NULL;
6481 }
6482 indices = PyMem_Malloc(n*sizeof(int));
6483 if (!indices) {
6484 PyErr_NoMemory();
6485 return NULL;
6486 }
6487 if (!mapping) for (i = 0; i < n; i++) indices[i] = s[i];
6488 else {
6489 for (i = 0; i < n; i++) {
6490 c = s[i];
6491 index = mapping[(int)c];
6492 if (index == MISSING_LETTER) {
6493 PyErr_SetString(PyExc_ValueError,
6494 "sequence contains letters not in the alphabet");
6495 PyMem_Free(indices);
6496 return NULL;
6497 }
6498 indices[i] = index;
6499 }
6500 }
6501 return indices;
6502 }
6503
6504 static int
convert_objects_to_ints(Py_buffer * view,PyObject * alphabet,PyObject * sequence)6505 convert_objects_to_ints(Py_buffer* view, PyObject* alphabet, PyObject* sequence)
6506 {
6507 Py_ssize_t i, j;
6508 Py_ssize_t n;
6509 Py_ssize_t m;
6510 int* indices = NULL;
6511 PyObject *obj1, *obj2;
6512 int equal;
6513
6514 view->buf = NULL;
6515 sequence = PySequence_Fast(sequence,
6516 "argument should support the sequence protocol");
6517 if (!sequence) return 0;
6518 if (!alphabet) {
6519 PyErr_SetString(PyExc_ValueError,
6520 "alphabet is None; cannot interpret sequence");
6521 goto exit;
6522 }
6523 alphabet = PySequence_Fast(alphabet, NULL); /* should never fail */
6524 n = PySequence_Size(sequence);
6525 m = PySequence_Size(alphabet);
6526 indices = PyMem_Malloc(n*sizeof(int));
6527 if (!indices) {
6528 PyErr_NoMemory();
6529 goto exit;
6530 }
6531 for (i = 0; i < n; i++) {
6532 obj1 = PySequence_Fast_GET_ITEM(sequence, i);
6533 for (j = 0; j < m; j++) {
6534 obj2 = PySequence_Fast_GET_ITEM(alphabet, j);
6535 equal = PyObject_RichCompareBool(obj1, obj2, Py_EQ);
6536 if (equal == 1) /* obj1 == obj2 */ {
6537 indices[i] = j;
6538 break;
6539 }
6540 else if (equal == -1) /* error */ {
6541 PyMem_Del(indices);
6542 goto exit;
6543 }
6544 /* else (equal == 0) continue; */ /* not equal */
6545 }
6546 if (j == m) {
6547 PyErr_SetString(PyExc_ValueError, "failed to find object in alphabet");
6548 goto exit;
6549 }
6550 }
6551 view->buf = indices;
6552 view->itemsize = 1;
6553 view->len = n;
6554 exit:
6555 Py_DECREF(sequence);
6556 Py_XDECREF(alphabet);
6557 if (view->buf) return 1;
6558 return 0;
6559 }
6560
6561 static int
sequence_converter(PyObject * argument,void * pointer)6562 sequence_converter(PyObject* argument, void* pointer)
6563 {
6564 Py_buffer* view = pointer;
6565 Py_ssize_t i;
6566 Py_ssize_t n;
6567 int index;
6568 int* indices;
6569 const int flag = PyBUF_FORMAT | PyBUF_C_CONTIGUOUS;
6570 Aligner* aligner;
6571 int* mapping;
6572
6573 if (argument == NULL) {
6574 if (view->obj) PyBuffer_Release(view);
6575 else {
6576 indices = view->buf;
6577 PyMem_Free(indices);
6578 }
6579 return 1;
6580 }
6581
6582 aligner = (Aligner*)view->obj;
6583 view->obj = NULL;
6584
6585 if (PyObject_GetBuffer(argument, view, flag) == 0) {
6586 if (view->ndim != 1) {
6587 PyErr_Format(PyExc_ValueError,
6588 "sequence has incorrect rank (%d expected 1)", view->ndim);
6589 return 0;
6590 }
6591 n = view->len / view->itemsize;
6592 if (n == 0) {
6593 PyErr_SetString(PyExc_ValueError, "sequence has zero length");
6594 return 0;
6595 }
6596 if (strcmp(view->format, "c") == 0 || strcmp(view->format, "B") == 0) {
6597 if (view->itemsize != sizeof(char)) {
6598 PyErr_Format(PyExc_ValueError,
6599 "sequence has unexpected item byte size "
6600 "(%ld, expected %ld)", view->itemsize, sizeof(char));
6601 return 0;
6602 }
6603 indices = convert_1bytes_to_ints(aligner->mapping, n, view->buf);
6604 if (!indices) return 0;
6605 PyBuffer_Release(view);
6606 view->itemsize = 1;
6607 view->len = n;
6608 view->buf = indices;
6609 return Py_CLEANUP_SUPPORTED;
6610 }
6611 if (strcmp(view->format, "i") == 0 || strcmp(view->format, "l") == 0) {
6612 if (view->itemsize != sizeof(int)) {
6613 PyErr_Format(PyExc_ValueError,
6614 "sequence has unexpected item byte size "
6615 "(%ld, expected %ld)", view->itemsize, sizeof(int));
6616 return 0;
6617 }
6618 indices = view->buf;
6619 if (aligner->substitution_matrix.obj) {
6620 const Py_ssize_t m = aligner->substitution_matrix.shape[0];
6621 for (i = 0; i < n; i++) {
6622 index = indices[i];
6623 if (index < 0) {
6624 PyErr_Format(PyExc_ValueError,
6625 "sequence item %zd is negative (%d)",
6626 i, index);
6627 return 0;
6628 }
6629 if (index >= m) {
6630 PyErr_Format(PyExc_ValueError,
6631 "sequence item %zd is out of bound"
6632 " (%d, should be < %zd)", i, index, m);
6633 return 0;
6634 }
6635 }
6636 }
6637 return Py_CLEANUP_SUPPORTED;
6638 }
6639 PyErr_Format(PyExc_ValueError,
6640 "sequence has incorrect data type '%s'", view->format);
6641 return 0;
6642 }
6643 PyErr_Clear(); /* To clear the exception raised by PyObject_GetBuffer */
6644 mapping = aligner->mapping;
6645 if (PyUnicode_Check(argument)) {
6646 if (PyUnicode_READY(argument) == -1) return 0;
6647 n = PyUnicode_GET_LENGTH(argument);
6648 switch (PyUnicode_KIND(argument)) {
6649 case PyUnicode_1BYTE_KIND: {
6650 Py_UCS1* s = PyUnicode_1BYTE_DATA(argument);
6651 indices = convert_1bytes_to_ints(mapping, n, (unsigned char*)s);
6652 break;
6653 }
6654 case PyUnicode_2BYTE_KIND: {
6655 Py_UCS2* s = PyUnicode_2BYTE_DATA(argument);
6656 indices = convert_2bytes_to_ints(mapping, n, s);
6657 break;
6658 }
6659 case PyUnicode_4BYTE_KIND: {
6660 Py_UCS4* s = PyUnicode_4BYTE_DATA(argument);
6661 indices = convert_4bytes_to_ints(mapping, n, s);
6662 break;
6663 }
6664 case PyUnicode_WCHAR_KIND:
6665 default:
6666 PyErr_SetString(PyExc_ValueError, "could not interpret unicode data");
6667 return 0;
6668 }
6669 if (!indices) return 0;
6670 view->buf = indices;
6671 view->itemsize = 1;
6672 view->len = n;
6673 return Py_CLEANUP_SUPPORTED;
6674 }
6675
6676 if (!mapping) {
6677 if (!convert_objects_to_ints(view, aligner->alphabet, argument)) return 0;
6678 return Py_CLEANUP_SUPPORTED;
6679 }
6680
6681 PyErr_SetString(PyExc_ValueError, "sequence has unexpected format");
6682 return 0;
6683 }
6684
6685 static int
strand_converter(PyObject * argument,void * pointer)6686 strand_converter(PyObject* argument, void* pointer)
6687 {
6688 if (!PyUnicode_Check(argument)) goto error;
6689 if (PyUnicode_READY(argument) == -1) return 0;
6690 if (PyUnicode_GET_LENGTH(argument) == 1) {
6691 const Py_UCS4 ch = PyUnicode_READ_CHAR(argument, 0);
6692 if (ch < 128) {
6693 const char c = ch;
6694 if (ch == '+' || ch == '-') {
6695 *((char*)pointer) = c;
6696 return 1;
6697 }
6698 }
6699 }
6700 error:
6701 PyErr_SetString(PyExc_ValueError, "strand must be '+' or '-'");
6702 return 0;
6703 }
6704
6705 static const char Aligner_score__doc__[] = "calculates the alignment score";
6706
6707 static PyObject*
Aligner_score(Aligner * self,PyObject * args,PyObject * keywords)6708 Aligner_score(Aligner* self, PyObject* args, PyObject* keywords)
6709 {
6710 const int* sA;
6711 const int* sB;
6712 Py_ssize_t nA;
6713 Py_ssize_t nB;
6714 Py_buffer bA = {0};
6715 Py_buffer bB = {0};
6716 const Mode mode = self->mode;
6717 const Algorithm algorithm = _get_algorithm(self);
6718 char strand = '+';
6719 PyObject* result = NULL;
6720 PyObject* substitution_matrix = self->substitution_matrix.obj;
6721
6722 static char *kwlist[] = {"sequenceA", "sequenceB", "strand", NULL};
6723
6724 bA.obj = (PyObject*)self;
6725 bB.obj = (PyObject*)self;
6726 if(!PyArg_ParseTupleAndKeywords(args, keywords, "O&O&O&", kwlist,
6727 sequence_converter, &bA,
6728 sequence_converter, &bB,
6729 strand_converter, &strand))
6730 return NULL;
6731
6732 sA = bA.buf;
6733 nA = bA.len / bA.itemsize;
6734 sB = bB.buf;
6735 nB = bB.len / bB.itemsize;
6736
6737 switch (algorithm) {
6738 case NeedlemanWunschSmithWaterman:
6739 switch (mode) {
6740 case Global:
6741 if (substitution_matrix)
6742 result = Aligner_needlemanwunsch_score_matrix(self, sA, nA, sB, nB, strand);
6743 else
6744 result = Aligner_needlemanwunsch_score_compare(self, sA, nA, sB, nB, strand);
6745 break;
6746 case Local:
6747 if (substitution_matrix)
6748 result = Aligner_smithwaterman_score_matrix(self, sA, nA, sB, nB);
6749 else
6750 result = Aligner_smithwaterman_score_compare(self, sA, nA, sB, nB);
6751 break;
6752 }
6753 break;
6754 case Gotoh:
6755 switch (mode) {
6756 case Global:
6757 if (substitution_matrix)
6758 result = Aligner_gotoh_global_score_matrix(self, sA, nA, sB, nB, strand);
6759 else
6760 result = Aligner_gotoh_global_score_compare(self, sA, nA, sB, nB, strand);
6761 break;
6762 case Local:
6763 if (substitution_matrix)
6764 result = Aligner_gotoh_local_score_matrix(self, sA, nA, sB, nB);
6765 else
6766 result = Aligner_gotoh_local_score_compare(self, sA, nA, sB, nB);
6767 break;
6768 }
6769 break;
6770 case WatermanSmithBeyer:
6771 switch (mode) {
6772 case Global:
6773 if (substitution_matrix)
6774 result = Aligner_watermansmithbeyer_global_score_matrix(self, sA, nA, sB, nB, strand);
6775 else
6776 result = Aligner_watermansmithbeyer_global_score_compare(self, sA, nA, sB, nB, strand);
6777 break;
6778 case Local:
6779 if (substitution_matrix)
6780 result = Aligner_watermansmithbeyer_local_score_matrix(self, sA, nA, sB, nB, strand);
6781 else
6782 result = Aligner_watermansmithbeyer_local_score_compare(self, sA, nA, sB, nB, strand);
6783 break;
6784 }
6785 break;
6786 case Unknown:
6787 default:
6788 PyErr_SetString(PyExc_RuntimeError, "unknown algorithm");
6789 break;
6790 }
6791
6792 sequence_converter(NULL, &bA);
6793 sequence_converter(NULL, &bB);
6794
6795 return result;
6796 }
6797
6798 static const char Aligner_align__doc__[] = "align two sequences";
6799
6800 static PyObject*
Aligner_align(Aligner * self,PyObject * args,PyObject * keywords)6801 Aligner_align(Aligner* self, PyObject* args, PyObject* keywords)
6802 {
6803 const int* sA;
6804 const int* sB;
6805 Py_ssize_t nA;
6806 Py_ssize_t nB;
6807 Py_buffer bA = {0};
6808 Py_buffer bB = {0};
6809 const Mode mode = self->mode;
6810 const Algorithm algorithm = _get_algorithm(self);
6811 char strand = '+';
6812 PyObject* result = NULL;
6813 PyObject* substitution_matrix = self->substitution_matrix.obj;
6814
6815 static char *kwlist[] = {"sequenceA", "sequenceB", "strand", NULL};
6816
6817 bA.obj = (PyObject*)self;
6818 bB.obj = (PyObject*)self;
6819 if(!PyArg_ParseTupleAndKeywords(args, keywords, "O&O&O&", kwlist,
6820 sequence_converter, &bA,
6821 sequence_converter, &bB,
6822 strand_converter, &strand))
6823 return NULL;
6824
6825 sA = bA.buf;
6826 nA = bA.len / bA.itemsize;
6827 sB = bB.buf;
6828 nB = bB.len / bB.itemsize;
6829
6830 switch (algorithm) {
6831 case NeedlemanWunschSmithWaterman:
6832 switch (mode) {
6833 case Global:
6834 if (substitution_matrix)
6835 result = Aligner_needlemanwunsch_align_matrix(self, sA, nA, sB, nB, strand);
6836 else
6837 result = Aligner_needlemanwunsch_align_compare(self, sA, nA, sB, nB, strand);
6838 break;
6839 case Local:
6840 if (substitution_matrix)
6841 result = Aligner_smithwaterman_align_matrix(self, sA, nA, sB, nB, strand);
6842 else
6843 result = Aligner_smithwaterman_align_compare(self, sA, nA, sB, nB, strand);
6844 break;
6845 }
6846 break;
6847 case Gotoh:
6848 switch (mode) {
6849 case Global:
6850 if (substitution_matrix)
6851 result = Aligner_gotoh_global_align_matrix(self, sA, nA, sB, nB, strand);
6852 else
6853 result = Aligner_gotoh_global_align_compare(self, sA, nA, sB, nB, strand);
6854 break;
6855 case Local:
6856 if (substitution_matrix)
6857 result = Aligner_gotoh_local_align_matrix(self, sA, nA, sB, nB, strand);
6858 else
6859 result = Aligner_gotoh_local_align_compare(self, sA, nA, sB, nB, strand);
6860 break;
6861 }
6862 break;
6863 case WatermanSmithBeyer:
6864 switch (mode) {
6865 case Global:
6866 if (substitution_matrix)
6867 result = Aligner_watermansmithbeyer_global_align_matrix(self, sA, nA, sB, nB, strand);
6868 else
6869 result = Aligner_watermansmithbeyer_global_align_compare(self, sA, nA, sB, nB, strand);
6870 break;
6871 case Local:
6872 if (substitution_matrix)
6873 result = Aligner_watermansmithbeyer_local_align_matrix(self, sA, nA, sB, nB, strand);
6874 else
6875 result = Aligner_watermansmithbeyer_local_align_compare(self, sA, nA, sB, nB, strand);
6876 break;
6877 }
6878 break;
6879 case Unknown:
6880 default:
6881 PyErr_SetString(PyExc_RuntimeError, "unknown algorithm");
6882 break;
6883 }
6884
6885 sequence_converter(NULL, &bA);
6886 sequence_converter(NULL, &bB);
6887
6888 return result;
6889 }
6890
6891 static char Aligner_doc[] =
6892 "Aligner.\n";
6893
6894 static PyMethodDef Aligner_methods[] = {
6895 {"score",
6896 (PyCFunction)Aligner_score,
6897 METH_VARARGS | METH_KEYWORDS,
6898 Aligner_score__doc__
6899 },
6900 {"align",
6901 (PyCFunction)Aligner_align,
6902 METH_VARARGS | METH_KEYWORDS,
6903 Aligner_align__doc__
6904 },
6905 {NULL} /* Sentinel */
6906 };
6907
6908 static PyTypeObject AlignerType = {
6909 PyVarObject_HEAD_INIT(NULL, 0)
6910 "_algorithms.PairwiseAligner", /* tp_name */
6911 sizeof(Aligner), /* tp_basicsize */
6912 0, /* tp_itemsize */
6913 (destructor)Aligner_dealloc, /* tp_dealloc */
6914 0, /* tp_print */
6915 0, /* tp_getattr */
6916 0, /* tp_setattr */
6917 0, /* tp_compare */
6918 (reprfunc)Aligner_repr, /* tp_repr */
6919 0, /* tp_as_number */
6920 0, /* tp_as_sequence */
6921 0, /* tp_as_mapping */
6922 0, /* tp_hash */
6923 0, /* tp_call */
6924 (reprfunc)Aligner_str, /* tp_str */
6925 0, /* tp_getattro */
6926 0, /* tp_setattro */
6927 0, /* tp_as_buffer */
6928 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
6929 Aligner_doc, /* tp_doc */
6930 0, /* tp_traverse */
6931 0, /* tp_clear */
6932 0, /* tp_richcompare */
6933 0, /* tp_weaklistoffset */
6934 0, /* tp_iter */
6935 0, /* tp_iternext */
6936 Aligner_methods, /* tp_methods */
6937 0, /* tp_members */
6938 Aligner_getset, /* tp_getset */
6939 0, /* tp_base */
6940 0, /* tp_dict */
6941 0, /* tp_descr_get */
6942 0, /* tp_descr_set */
6943 0, /* tp_dictoffset */
6944 (initproc)Aligner_init, /* tp_init */
6945 };
6946
6947
6948 /* Module definition */
6949
6950 static char _aligners__doc__[] =
6951 "C extension module implementing pairwise alignment algorithms";
6952
6953 static struct PyModuleDef moduledef = {
6954 PyModuleDef_HEAD_INIT,
6955 "_aligners",
6956 _aligners__doc__,
6957 -1,
6958 NULL,
6959 NULL,
6960 NULL,
6961 NULL,
6962 NULL
6963 };
6964
6965 PyObject *
PyInit__aligners(void)6966 PyInit__aligners(void)
6967 {
6968 PyObject* module;
6969 AlignerType.tp_new = PyType_GenericNew;
6970
6971 if (PyType_Ready(&AlignerType) < 0 || PyType_Ready(&PathGenerator_Type) < 0)
6972 return NULL;
6973
6974 module = PyModule_Create(&moduledef);
6975 if (!module) return NULL;
6976
6977 Py_INCREF(&AlignerType);
6978 /* Reference to AlignerType will be stolen by PyModule_AddObject
6979 * only if it is successful. */
6980 if (PyModule_AddObject(module,
6981 "PairwiseAligner", (PyObject*) &AlignerType) < 0) {
6982 Py_DECREF(&AlignerType);
6983 Py_DECREF(module);
6984 return NULL;
6985 }
6986
6987 return module;
6988 }
6989