jpayne@69: /* Copyright 2002 by Jeffrey Chang. jpayne@69: * Copyright 2016, 2019 by Markus Piotrowski. jpayne@69: * All rights reserved. jpayne@69: * jpayne@69: * This file is part of the Biopython distribution and governed by your jpayne@69: * choice of the "Biopython License Agreement" or the "BSD 3-Clause License". jpayne@69: * Please see the LICENSE file that should have been included as part of this jpayne@69: * package. jpayne@69: * jpayne@69: * cpairwise2module.c jpayne@69: * Created 30 Sep 2001 jpayne@69: * jpayne@69: * Optimized C routines that complement pairwise2.py. jpayne@69: */ jpayne@69: jpayne@69: #include "Python.h" jpayne@69: jpayne@69: jpayne@69: #define _PRECISION 1000 jpayne@69: #define rint(x) (int)((x)*_PRECISION+0.5) jpayne@69: jpayne@69: /* Functions in this module. */ jpayne@69: jpayne@69: static double calc_affine_penalty(int length, double open, double extend, jpayne@69: int penalize_extend_when_opening) jpayne@69: { jpayne@69: double penalty; jpayne@69: jpayne@69: if(length <= 0) jpayne@69: return 0.0; jpayne@69: penalty = open + extend * length; jpayne@69: if(!penalize_extend_when_opening) jpayne@69: penalty -= extend; jpayne@69: return penalty; jpayne@69: } jpayne@69: jpayne@69: static double _get_match_score(PyObject *py_sequenceA, PyObject *py_sequenceB, jpayne@69: PyObject *py_match_fn, int i, int j, jpayne@69: char *sequenceA, char *sequenceB, jpayne@69: int use_sequence_cstring, jpayne@69: double match, double mismatch, jpayne@69: int use_match_mismatch_scores) jpayne@69: { jpayne@69: PyObject *py_A=NULL, *py_B=NULL; jpayne@69: PyObject *py_arglist=NULL, *py_result=NULL; jpayne@69: double score = 0; jpayne@69: jpayne@69: if(use_sequence_cstring && use_match_mismatch_scores) { jpayne@69: score = (sequenceA[i] == sequenceB[j]) ? match : mismatch; jpayne@69: return score; jpayne@69: } jpayne@69: /* Calculate the match score. */ jpayne@69: if(!(py_A = PySequence_GetItem(py_sequenceA, i))) jpayne@69: goto _get_match_score_cleanup; jpayne@69: if(!(py_B = PySequence_GetItem(py_sequenceB, j))) jpayne@69: goto _get_match_score_cleanup; jpayne@69: if(!(py_arglist = Py_BuildValue("(OO)", py_A, py_B))) jpayne@69: goto _get_match_score_cleanup; jpayne@69: jpayne@69: if(!(py_result = PyEval_CallObject(py_match_fn, py_arglist))) jpayne@69: goto _get_match_score_cleanup; jpayne@69: score = PyFloat_AsDouble(py_result); jpayne@69: jpayne@69: _get_match_score_cleanup: jpayne@69: if(py_A) { jpayne@69: Py_DECREF(py_A); jpayne@69: } jpayne@69: if(py_B) { jpayne@69: Py_DECREF(py_B); jpayne@69: } jpayne@69: if(py_arglist) { jpayne@69: Py_DECREF(py_arglist); jpayne@69: } jpayne@69: if(py_result) { jpayne@69: Py_DECREF(py_result); jpayne@69: } jpayne@69: return score; jpayne@69: } jpayne@69: jpayne@69: #if PY_MAJOR_VERSION >= 3 jpayne@69: static PyObject* _create_bytes_object(PyObject* o) jpayne@69: { jpayne@69: PyObject* b; jpayne@69: if (PyBytes_Check(o)) { jpayne@69: return o; jpayne@69: } jpayne@69: if (!PyUnicode_Check(o)) { jpayne@69: return NULL; jpayne@69: } jpayne@69: b = PyUnicode_AsASCIIString(o); jpayne@69: if (!b) { jpayne@69: PyErr_Clear(); jpayne@69: return NULL; jpayne@69: } jpayne@69: return b; jpayne@69: } jpayne@69: #endif jpayne@69: jpayne@69: /* This function is a more-or-less straightforward port of the jpayne@69: * equivalent function in pairwise2. Please see there for algorithm jpayne@69: * documentation. jpayne@69: */ jpayne@69: static PyObject *cpairwise2__make_score_matrix_fast(PyObject *self, jpayne@69: PyObject *args) jpayne@69: { jpayne@69: int i; jpayne@69: int row, col; jpayne@69: PyObject *py_sequenceA, *py_sequenceB, *py_match_fn; jpayne@69: #if PY_MAJOR_VERSION >= 3 jpayne@69: PyObject *py_bytesA, *py_bytesB; jpayne@69: #endif jpayne@69: char *sequenceA=NULL, *sequenceB=NULL; jpayne@69: int use_sequence_cstring; jpayne@69: double open_A, extend_A, open_B, extend_B; jpayne@69: int penalize_extend_when_opening, penalize_end_gaps_A, penalize_end_gaps_B; jpayne@69: int align_globally, score_only; jpayne@69: jpayne@69: PyObject *py_match=NULL, *py_mismatch=NULL; jpayne@69: double first_A_gap, first_B_gap; jpayne@69: double match, mismatch; jpayne@69: double score; jpayne@69: double best_score = 0; jpayne@69: double local_max_score = 0; jpayne@69: int use_match_mismatch_scores; jpayne@69: int lenA, lenB; jpayne@69: double *score_matrix = NULL; jpayne@69: unsigned char *trace_matrix = NULL; jpayne@69: PyObject *py_score_matrix=NULL, *py_trace_matrix=NULL; jpayne@69: jpayne@69: double *col_cache_score = NULL; jpayne@69: PyObject *py_retval = NULL; jpayne@69: jpayne@69: if(!PyArg_ParseTuple(args, "OOOddddi(ii)ii", &py_sequenceA, &py_sequenceB, jpayne@69: &py_match_fn, &open_A, &extend_A, &open_B, &extend_B, jpayne@69: &penalize_extend_when_opening, jpayne@69: &penalize_end_gaps_A, &penalize_end_gaps_B, jpayne@69: &align_globally, &score_only)) jpayne@69: return NULL; jpayne@69: if(!PySequence_Check(py_sequenceA) || !PySequence_Check(py_sequenceB)) { jpayne@69: PyErr_SetString(PyExc_TypeError, jpayne@69: "py_sequenceA and py_sequenceB should be sequences."); jpayne@69: return NULL; jpayne@69: } jpayne@69: jpayne@69: /* Optimize for the common case. Check to see if py_sequenceA and jpayne@69: py_sequenceB are strings. If they are, use the c string jpayne@69: representation. */ jpayne@69: #if PY_MAJOR_VERSION < 3 jpayne@69: use_sequence_cstring = 0; jpayne@69: if(PyString_Check(py_sequenceA) && PyString_Check(py_sequenceB)) { jpayne@69: sequenceA = PyString_AS_STRING(py_sequenceA); jpayne@69: sequenceB = PyString_AS_STRING(py_sequenceB); jpayne@69: use_sequence_cstring = 1; jpayne@69: } jpayne@69: #else jpayne@69: py_bytesA = _create_bytes_object(py_sequenceA); jpayne@69: py_bytesB = _create_bytes_object(py_sequenceB); jpayne@69: if (py_bytesA && py_bytesB) { jpayne@69: sequenceA = PyBytes_AS_STRING(py_bytesA); jpayne@69: sequenceB = PyBytes_AS_STRING(py_bytesB); jpayne@69: use_sequence_cstring = 1; jpayne@69: } jpayne@69: else { jpayne@69: Py_XDECREF(py_bytesA); jpayne@69: Py_XDECREF(py_bytesB); jpayne@69: use_sequence_cstring = 0; jpayne@69: } jpayne@69: #endif jpayne@69: jpayne@69: if(!PyCallable_Check(py_match_fn)) { jpayne@69: PyErr_SetString(PyExc_TypeError, "py_match_fn must be callable."); jpayne@69: return NULL; jpayne@69: } jpayne@69: /* Optimize for the common case. Check to see if py_match_fn is jpayne@69: an identity_match. If so, pull out the match and mismatch jpayne@69: member variables and calculate the scores myself. */ jpayne@69: match = mismatch = 0; jpayne@69: use_match_mismatch_scores = 0; jpayne@69: if(!(py_match = PyObject_GetAttrString(py_match_fn, "match"))) jpayne@69: goto cleanup_after_py_match_fn; jpayne@69: match = PyFloat_AsDouble(py_match); jpayne@69: if(match==-1.0 && PyErr_Occurred()) jpayne@69: goto cleanup_after_py_match_fn; jpayne@69: if(!(py_mismatch = PyObject_GetAttrString(py_match_fn, "mismatch"))) jpayne@69: goto cleanup_after_py_match_fn; jpayne@69: mismatch = PyFloat_AsDouble(py_mismatch); jpayne@69: if(mismatch==-1.0 && PyErr_Occurred()) jpayne@69: goto cleanup_after_py_match_fn; jpayne@69: use_match_mismatch_scores = 1; jpayne@69: jpayne@69: cleanup_after_py_match_fn: jpayne@69: if(PyErr_Occurred()) jpayne@69: PyErr_Clear(); jpayne@69: if(py_match) { jpayne@69: Py_DECREF(py_match); jpayne@69: } jpayne@69: if(py_mismatch) { jpayne@69: Py_DECREF(py_mismatch); jpayne@69: } jpayne@69: /* Cache some commonly used gap penalties */ jpayne@69: first_A_gap = calc_affine_penalty(1, open_A, extend_A, jpayne@69: penalize_extend_when_opening); jpayne@69: first_B_gap = calc_affine_penalty(1, open_B, extend_B, jpayne@69: penalize_extend_when_opening); jpayne@69: jpayne@69: /* Allocate matrices for storing the results and initialize first row and col. */ jpayne@69: lenA = PySequence_Length(py_sequenceA); jpayne@69: lenB = PySequence_Length(py_sequenceB); jpayne@69: score_matrix = malloc((lenA+1)*(lenB+1)*sizeof(*score_matrix)); jpayne@69: if(!score_matrix) { jpayne@69: PyErr_SetString(PyExc_MemoryError, "Out of memory"); jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: } jpayne@69: for(i=0; i<(lenB+1); i++) jpayne@69: score_matrix[i] = 0; jpayne@69: for(i=0; i<(lenA+1)*(lenB+1); i += (lenB+1)) jpayne@69: score_matrix[i] = 0; jpayne@69: /* If we only want the score, we don't need the trace matrix. */ jpayne@69: if (!score_only){ jpayne@69: trace_matrix = malloc((lenA+1)*(lenB+1)*sizeof(*trace_matrix)); jpayne@69: if(!trace_matrix) { jpayne@69: PyErr_SetString(PyExc_MemoryError, "Out of memory"); jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: } jpayne@69: for(i=0; i<(lenB+1); i++) jpayne@69: trace_matrix[i] = 0; jpayne@69: for(i=0; i<(lenA+1)*(lenB+1); i += (lenB+1)) jpayne@69: trace_matrix[i] = 0; jpayne@69: } jpayne@69: else jpayne@69: trace_matrix = malloc(1); jpayne@69: jpayne@69: /* Initialize the first row and col of the score matrix. */ jpayne@69: for(i=0; i<=lenA; i++) { jpayne@69: if(penalize_end_gaps_B) jpayne@69: score = calc_affine_penalty(i, open_B, extend_B, jpayne@69: penalize_extend_when_opening); jpayne@69: else jpayne@69: score = 0; jpayne@69: score_matrix[i*(lenB+1)] = score; jpayne@69: } jpayne@69: for(i=0; i<=lenB; i++) { jpayne@69: if(penalize_end_gaps_A) jpayne@69: score = calc_affine_penalty(i, open_A, extend_A, jpayne@69: penalize_extend_when_opening); jpayne@69: else jpayne@69: score = 0; jpayne@69: score_matrix[i] = score; jpayne@69: } jpayne@69: jpayne@69: /* Now initialize the col cache. */ jpayne@69: col_cache_score = malloc((lenB+1)*sizeof(*col_cache_score)); jpayne@69: memset((void *)col_cache_score, 0, (lenB+1)*sizeof(*col_cache_score)); jpayne@69: for(i=0; i<=lenB; i++) { jpayne@69: col_cache_score[i] = calc_affine_penalty(i, (2*open_B), extend_B, jpayne@69: penalize_extend_when_opening); jpayne@69: } jpayne@69: jpayne@69: /* Fill in the score matrix. The row cache is calculated on the fly.*/ jpayne@69: for(row=1; row<=lenA; row++) { jpayne@69: double row_cache_score = calc_affine_penalty(row, (2*open_A), extend_A, jpayne@69: penalize_extend_when_opening); jpayne@69: for(col=1; col<=lenB; col++) { jpayne@69: double match_score, nogap_score; jpayne@69: double row_open, row_extend, col_open, col_extend; jpayne@69: int best_score_rint, row_score_rint, col_score_rint; jpayne@69: unsigned char row_trace_score, col_trace_score, trace_score; jpayne@69: jpayne@69: /* Calculate the best score. */ jpayne@69: match_score = _get_match_score(py_sequenceA, py_sequenceB, jpayne@69: py_match_fn, row-1, col-1, jpayne@69: sequenceA, sequenceB, jpayne@69: use_sequence_cstring, jpayne@69: match, mismatch, jpayne@69: use_match_mismatch_scores); jpayne@69: if(match_score==-1.0 && PyErr_Occurred()) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: nogap_score = score_matrix[(row-1)*(lenB+1)+col-1] + match_score; jpayne@69: jpayne@69: if (!penalize_end_gaps_A && row==lenA) { jpayne@69: row_open = score_matrix[(row)*(lenB+1)+col-1]; jpayne@69: row_extend = row_cache_score; jpayne@69: } jpayne@69: else { jpayne@69: row_open = score_matrix[(row)*(lenB+1)+col-1] + first_A_gap; jpayne@69: row_extend = row_cache_score + extend_A; jpayne@69: } jpayne@69: row_cache_score = (row_open > row_extend) ? row_open : row_extend; jpayne@69: jpayne@69: if (!penalize_end_gaps_B && col==lenB){ jpayne@69: col_open = score_matrix[(row-1)*(lenB+1)+col]; jpayne@69: col_extend = col_cache_score[col]; jpayne@69: } jpayne@69: else { jpayne@69: col_open = score_matrix[(row-1)*(lenB+1)+col] + first_B_gap; jpayne@69: col_extend = col_cache_score[col] + extend_B; jpayne@69: } jpayne@69: col_cache_score[col] = (col_open > col_extend) ? col_open : col_extend; jpayne@69: jpayne@69: best_score = (row_cache_score > col_cache_score[col]) ? row_cache_score : col_cache_score[col]; jpayne@69: if(nogap_score > best_score) jpayne@69: best_score = nogap_score; jpayne@69: jpayne@69: if (best_score > local_max_score) jpayne@69: local_max_score = best_score; jpayne@69: jpayne@69: if(!align_globally && best_score < 0) jpayne@69: score_matrix[row*(lenB+1)+col] = 0; jpayne@69: else jpayne@69: score_matrix[row*(lenB+1)+col] = best_score; jpayne@69: jpayne@69: if (!score_only) { jpayne@69: row_score_rint = rint(row_cache_score); jpayne@69: col_score_rint = rint(col_cache_score[col]); jpayne@69: row_trace_score = 0; jpayne@69: col_trace_score = 0; jpayne@69: if (rint(row_open) == row_score_rint) jpayne@69: row_trace_score = row_trace_score|1; jpayne@69: if (rint(row_extend) == row_score_rint) jpayne@69: row_trace_score = row_trace_score|8; jpayne@69: if (rint(col_open) == col_score_rint) jpayne@69: col_trace_score = col_trace_score|4; jpayne@69: if (rint(col_extend) == col_score_rint) jpayne@69: col_trace_score = col_trace_score|16; jpayne@69: jpayne@69: trace_score = 0; jpayne@69: best_score_rint = rint(best_score); jpayne@69: if (rint(nogap_score) == best_score_rint) jpayne@69: trace_score = trace_score|2; jpayne@69: if (row_score_rint == best_score_rint) jpayne@69: trace_score += row_trace_score; jpayne@69: if (col_score_rint == best_score_rint) jpayne@69: trace_score += col_trace_score; jpayne@69: trace_matrix[row*(lenB+1)+col] = trace_score; jpayne@69: } jpayne@69: } jpayne@69: } jpayne@69: jpayne@69: if (!align_globally) jpayne@69: best_score = local_max_score; jpayne@69: jpayne@69: /* Save the score and traceback matrices into real python objects. */ jpayne@69: if(!score_only) { jpayne@69: if(!(py_score_matrix = PyList_New(lenA+1))) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: if(!(py_trace_matrix = PyList_New(lenA+1))) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: jpayne@69: for(row=0; row<=lenA; row++) { jpayne@69: PyObject *py_score_row, *py_trace_row; jpayne@69: if(!(py_score_row = PyList_New(lenB+1))) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: PyList_SET_ITEM(py_score_matrix, row, py_score_row); jpayne@69: if(!(py_trace_row = PyList_New(lenB+1))) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: PyList_SET_ITEM(py_trace_matrix, row, py_trace_row); jpayne@69: jpayne@69: for(col=0; col<=lenB; col++) { jpayne@69: PyObject *py_score, *py_trace; jpayne@69: int offset = row*(lenB+1) + col; jpayne@69: jpayne@69: /* Set py_score_matrix[row][col] to the score. */ jpayne@69: if(!(py_score = PyFloat_FromDouble(score_matrix[offset]))) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: PyList_SET_ITEM(py_score_row, col, py_score); jpayne@69: jpayne@69: /* Set py_trace_matrix[row][col] to a list of indexes. On jpayne@69: the edges of the matrix (row or column is 0), the jpayne@69: matrix should be [None]. */ jpayne@69: if(!row || !col) { jpayne@69: if(!(py_trace = Py_BuildValue("B", 1))) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: Py_INCREF(Py_None); jpayne@69: PyList_SET_ITEM(py_trace_row, col, Py_None); jpayne@69: } jpayne@69: else { jpayne@69: if(!(py_trace = Py_BuildValue("B", trace_matrix[offset]))) jpayne@69: goto _cleanup_make_score_matrix_fast; jpayne@69: PyList_SET_ITEM(py_trace_row, col, py_trace); jpayne@69: jpayne@69: } jpayne@69: } jpayne@69: } jpayne@69: } jpayne@69: else { jpayne@69: py_score_matrix = PyList_New(1); jpayne@69: py_trace_matrix = PyList_New(1); jpayne@69: } jpayne@69: py_retval = Py_BuildValue("(OOd)", py_score_matrix, py_trace_matrix, best_score); jpayne@69: jpayne@69: _cleanup_make_score_matrix_fast: jpayne@69: if(score_matrix) jpayne@69: free(score_matrix); jpayne@69: if(trace_matrix) jpayne@69: free(trace_matrix); jpayne@69: if(col_cache_score) jpayne@69: free(col_cache_score); jpayne@69: if(py_score_matrix){ jpayne@69: Py_DECREF(py_score_matrix); jpayne@69: } jpayne@69: if(py_trace_matrix){ jpayne@69: Py_DECREF(py_trace_matrix); jpayne@69: } jpayne@69: jpayne@69: #if PY_MAJOR_VERSION >= 3 jpayne@69: if (py_bytesA != NULL && py_bytesA != py_sequenceA) Py_DECREF(py_bytesA); jpayne@69: if (py_bytesB != NULL && py_bytesB != py_sequenceB) Py_DECREF(py_bytesB); jpayne@69: #endif jpayne@69: jpayne@69: return py_retval; jpayne@69: } jpayne@69: jpayne@69: static PyObject *cpairwise2_rint(PyObject *self, PyObject *args, jpayne@69: PyObject *keywds) jpayne@69: { jpayne@69: double x; jpayne@69: int precision = _PRECISION; jpayne@69: int rint_x; jpayne@69: jpayne@69: static char *kwlist[] = {"x", "precision", NULL}; jpayne@69: jpayne@69: if(!PyArg_ParseTupleAndKeywords(args, keywds, "d|l", kwlist, jpayne@69: &x, &precision)) jpayne@69: return NULL; jpayne@69: rint_x = (int)(x * precision + 0.5); jpayne@69: #if PY_MAJOR_VERSION >= 3 jpayne@69: return PyLong_FromLong((long)rint_x); jpayne@69: #else jpayne@69: return PyInt_FromLong((long)rint_x); jpayne@69: #endif jpayne@69: } jpayne@69: jpayne@69: /* Module definition stuff */ jpayne@69: jpayne@69: static PyMethodDef cpairwise2Methods[] = { jpayne@69: {"_make_score_matrix_fast", jpayne@69: (PyCFunction)cpairwise2__make_score_matrix_fast, METH_VARARGS, ""}, jpayne@69: {"rint", (PyCFunction)cpairwise2_rint, METH_VARARGS|METH_KEYWORDS, ""}, jpayne@69: {NULL, NULL, 0, NULL} jpayne@69: }; jpayne@69: jpayne@69: static char cpairwise2__doc__[] = jpayne@69: "Optimized C routines that complement pairwise2.py. These are called from within pairwise2.py.\n\ jpayne@69: \n\ jpayne@69: "; jpayne@69: jpayne@69: #if PY_MAJOR_VERSION >= 3 jpayne@69: jpayne@69: static struct PyModuleDef moduledef = { jpayne@69: PyModuleDef_HEAD_INIT, jpayne@69: "cpairwise2", jpayne@69: cpairwise2__doc__, jpayne@69: -1, jpayne@69: cpairwise2Methods, jpayne@69: NULL, jpayne@69: NULL, jpayne@69: NULL, jpayne@69: NULL jpayne@69: }; jpayne@69: jpayne@69: PyObject * jpayne@69: PyInit_cpairwise2(void) jpayne@69: jpayne@69: #else jpayne@69: jpayne@69: void jpayne@69: /* for Windows: _declspec(dllexport) initcpairwise2(void) */ jpayne@69: initcpairwise2(void) jpayne@69: #endif jpayne@69: jpayne@69: { jpayne@69: #if PY_MAJOR_VERSION >= 3 jpayne@69: PyObject* module = PyModule_Create(&moduledef); jpayne@69: if (module==NULL) return NULL; jpayne@69: return module; jpayne@69: #else jpayne@69: (void) Py_InitModule3("cpairwise2", cpairwise2Methods, cpairwise2__doc__); jpayne@69: #endif jpayne@69: }