//
//  GlobalAlignment.m
//  Alignment
//
//  Created by Alexander Griekspoor on 25-12-04.
//  Copyright 2004 Mekentosj.com. All rights reserved.
//

#import "GlobalAlignment.h"
#import "NSString_Extension.h"

// Directions
typedef enum {
    kNone = 0,
	kDiagonal,
    kLeft,
    kUp
} Pointers;

// Gaps
typedef enum {
	kNoGap = 0,
    kGapOpened
} Gaps;



int ** twodim_alloc(int length_x, int length_y);


@implementation GlobalAlignment



+ (NSDictionary *)alignSequence: (NSString *)sequence1 
				   withSequence: (NSString *)sequence2 
						  match: (int)match
					   mismatch: (int)mismatch
							gap: (int)gap {
	
	//  NEEDLEMAN - WUNSCH GLOBAL ALIGNMENT (BASIC)
	//
	//	    W   O   R   D   1  (A, M, i)
	//	
	//	W   0   0   0   0   0
	//	
	//	O   0   1	
	//	
	//	R   0
	// 	
	//	D   0
	//	
	//	2   0
	//	
	//	(B, N, j)
	
	int i, j, k;		// general counters
	int M, N;			// length of sequences
	int diag, left, up;	// temp scores
	char  *A,  *B;		// c strings equivalent of NSStrings
	char  *a,  *b;		// resulting strings
	int **score;		// score matrix
	int **pointer;		// pointer matrix
	
	// PREPARATION
	// size of sequences
	M = [sequence1 cStringLength];
	N = [sequence2 cStringLength];

	// allocate space for A and read in 
	A = ( char * ) malloc( (M + 1) * sizeof(char));
	[sequence1 getCString: A];
	
	// allocate space for B and read in 
	B = ( char * ) malloc( (N + 1) * sizeof(char));
	[sequence2 getCString: B];
	
	// INITIALIZATION
	// declare matrices to holds scores and pointers
	// first = rows, second = column
	int seqlen1 = M+1;
	int seqlen2 = N+1;
	score = twodim_alloc(seqlen1, seqlen2); 
	pointer = twodim_alloc(seqlen1, seqlen2); 
	
	// check
	if(score == NULL || pointer == NULL){
		NSLog(@"Alignment: could not allocate memory.");
		return nil;
	}
	
	// zero origin
	score[0][0] = 0;
	pointer[0][0] = kNone;
	
	// top row + left column
	for ( i = 1; i < M + 1 ; i++ ){
		score[i][0] = gap * i;
		pointer[i][0] = kLeft;
	}
	for ( j = 1; j < N + 1; j++ ){
		  score[0][j] = gap * j;
		pointer[0][j] = kUp;
	}

	// FILL
	// main loop, for every row
	for ( i = 1; i < M + 1 ; i++ ){
		// for every column
		for ( j = 1; j < N + 1 ; j++ ){

			// calculate diagonal score
			if(A[i-1] == B[j-1]) diag = score[i-1][j-1] + match;
			else diag = score[i-1][j-1] + mismatch;

			// calculate gap scores
			left = score[i - 1][j] + gap;
			  up = score[i][j - 1] + gap;
			
			// choose best scores
			if(diag >= up){
				if(diag >= left) {
					  score[i][j] = diag;
					pointer[i][j] = kDiagonal;
				}
				else {
					  score[i][j] = left;
					pointer[i][j] = kLeft;
				}
			} else {
				if(up >= left) {
					  score[i][j] = up;
					pointer[i][j] = kUp;
				}
				else {
					  score[i][j] = left;
					pointer[i][j] = kLeft;
				}	
			}
			//NSLog(@"i=%d, j=%d, score: %d, pointer: %d", i, j, score[i][j], pointer[i][j]);
		}
	}
	

	// TRACE BACK
	i = M;
	j = N;
	k = 0;
	a = ( char * ) malloc( (M + N + 1) * sizeof(char));
	b = ( char * ) malloc( (N + N + 1) * sizeof(char));

	while(1){
		// escape when origin is reached
		if(pointer[i][j] == kNone) break;
		
		switch(pointer[i][j]){
			case kDiagonal :
				a[k] = A[i - 1];
				b[k] = B[j - 1];
				i--;
				j--;
				k++;
				break;
				
			case kLeft :
				a[k] = A[i - 1];
				b[k] = '-';
				i--;
				k++;
				break;
				
			case kUp :
				a[k] = '-';
				b[k] = B[j - 1];
				j--;
				k++;
				break;
		}
	}

	
	// Prepare for export
	NSString *s1 = [[NSString stringWithCString: a length: k] convertToReverse];
	NSString *s2 = [[NSString stringWithCString: b length: k] convertToReverse];
	NSDictionary *dict = [NSDictionary dictionaryWithObjectsAndKeys: s1, @"Sequence 1", s2, @"Sequence 2", NULL];
	
	for(i=k-1;i>=0;i--) printf("%c",a[i]);
	printf("\n");
	for(j=k-1;j>=0;j--) printf("%c",b[j]);
	printf("\n");
	
	return dict;
}




@end



int ** twodim_alloc(int length_x, int length_y) 
{
	int **matrix;
	unsigned y;
	
	if (!(matrix = (int**) malloc(sizeof(int*) * (size_t) length_y)))	return NULL;				//Exit("memory allocation failed.", "", 0);
	if (!(*matrix = (int*) malloc(sizeof(int) * (size_t) (length_x * length_y))))	return NULL;	//Exit("memory allocation failed.", "", 0);
	for (y=1; y<length_y; y++)	matrix[y] = matrix[y-1] + length_x;
	return(matrix);
}