biohazard-2.0: cbits/myers_align.c

#include "myers_align.h"

#include <limits.h>
#include <stdlib.h>
#include <string.h>

// [*blech*, this looks and feels like FORTRAN.]
unsigned myers_diff(
        const char *seq_a, int len_a, enum myers_align_mode mode, 
        const char* seq_b, int len_b, int maxd,
        char *bt_a, char *bt_b ) 
{
	// int len_a = strlen( seq_a ), len_b = strlen( seq_b ) ;
	if( maxd > len_a + len_b ) maxd = len_a + len_b ;

	// in vee[d][k], d runs from 0 to maxd; k runs from -d to +d
	int **vee = calloc( maxd, sizeof(int*) ) ;

	int d, dd, k, x, y, r = UINT_MAX ;
	int *v_d_1 = 0, *v_d = 0 ; 															// "array slice" vee[.][d-1]
	for( d = 0 ; d != maxd ; ++d, v_d_1 = v_d )									// D-paths in order of increasing D
	{
		v_d = d + (vee[d] = malloc( (2 * d + 1) * sizeof( int ) )) ; 		// "array slice" vee[.][d]

		for( k = max(-d,-len_a) ; k <= min(d,len_b) ; ++k ) 					// diagonals
		{
			if( d == 0 )         x = 0 ;
			else if(d==1&&k==0)  x =                       v_d_1[ k ]+1 ;
			else if( k == -d   ) x =                                     v_d_1[ k+1 ] ;
			else if( k ==  d   ) x =       v_d_1[ k-1 ]+1 ;									// argh, need to check for d first, b/c -d+2 could be equal to d
			else if( k == -d+1 ) x = max(                  v_d_1[ k ]+1, v_d_1[ k+1 ] ) ;
			else if( k ==  d-1 ) x = max(  v_d_1[ k-1 ]+1, v_d_1[ k ]+1 ) ;
			else                 x = max3( v_d_1[ k-1 ]+1, v_d_1[ k ]+1, v_d_1[ k+1 ] ) ;

			y = x-k ;
			while( x < len_b && y < len_a && match( seq_b[x], seq_a[y] ) ) ++x, ++y ;
			v_d[ k ] = x ;

			if(
                    bt_a && bt_b &&
					(mode == myers_align_is_prefix || y == len_a) &&
					(mode == myers_align_has_prefix || x == len_b) )
			{
				char *out_a = bt_a + len_a + d +2 ;
				char *out_b = bt_b + len_b + d +2 ;
				*--out_a = 0 ;
				*--out_b = 0 ;
				for( dd = d ; dd != 0 ; )
				{
					if( k != -dd && k != dd && x == vee[ dd-1 ][ k + dd-1 ]+1 )
					{
						--dd ;
						--x ;
						--y ;
						*--out_b = seq_b[x] ;
						*--out_a = seq_a[y] ;
					}
					else if( k > -dd+1 && x == vee[ dd-1 ][ k-1 + dd-1 ]+1 )
					{
						--x ;
						--k ;
						--dd ;
						*--out_b = seq_b[x] ;
						*--out_a = '-' ;
					}
					else if( k < dd-1 && x == vee[ dd-1 ][ k+1 + dd-1 ] )
					{
						++k ;
						--y ;
						--dd ;
						*--out_b = '-' ;
						*--out_a = seq_a[y] ;
					}
					else // this better had been a match...
					{
						--x ;
						--y ;
						*--out_b = seq_b[x] ;
						*--out_a = seq_a[y] ;
					}
				}
				while( x > 0 )
				{
					--x ;
					*--out_b = seq_b[x] ;
					*--out_a = seq_a[x] ;
				}
				memmove( bt_a, out_a, bt_a + len_a + d + 2 - out_a ) ;
				memmove( bt_b, out_b, bt_b + len_b + d + 2 - out_b ) ;
				r = d ;
				goto cleanup ;
			}
		}
	}

cleanup:
	for( dd = maxd ; dd != 0 ; --dd )
		free( vee[dd-1] ) ;
	free( vee ) ;
	return r ;
}