#include <iostream>
#include <cstring> 
using namespace std;

int Max(int x, int y){
	if(x > y) return x;
	else return y;
}

int main () {
	string seq1;
	cout << "Please input the first DNA sequence: ";
	getline (cin, seq1);

	string seq2;
	cout << "Please input the second DNA sequence: ";
	getline (cin, seq2);

	char gl[3];
	cout << "Global or local alignment? (g/l) ";
	cin >> gl;

	int match;
	cout << "Please input the match score: ";
	cin >> match;

	int mismatch;
	cout << "Please input the mismatch score: ";
	cin >> mismatch;

	int indel;
	cout << "Please input the gap penalty: ";
	cin >> indel;

	bool type = (strcmp(gl, "g")==0 || strcmp(gl,"G")==0); // to determine global or local alignment
	int L1, L2;
	L1 = seq1.size(); 
	L2 = seq2.size();
	int V[300][300], Max_value = 0, Max_r = L1, Max_s = L2;
	char *align1, *align2; 
	align1 = new char[L1+L2];
	align2 = new char[L1+L2];
	
	// Boundary conditions
	
	if(type){// for global alignment 
		for(int j = L2;j>=0;j--) V[L1][j] = (L2-j)*indel;
		for(int i = L1;i>=0;i--) V[i][L2] = (L1-i)*indel;
	}else{// for local alignment 
		for(int j = L2;j>=0;j--) V[L1][j] = 0;
		for(int i = L1;i>=0;i--) V[i][L2] = 0;
	}
	
	// Recurrence relation
	
	for(int i = L1-1; i>=0;i--){
		for(int j = L2-1;j>=0;j--){
			V[i][j] = Max(V[i][j+1]+indel, V[i+1][j]+indel);
			if(seq1[i] == seq2[j]) V[i][j] = Max(V[i][j], V[i+1][j+1]+match);
			else V[i][j] = Max(V[i][j], V[i+1][j+1]+mismatch);
			if(type == false){
				V[i][j] = Max(V[i][j],0); // extra condition for local alignment
				if(V[i][j] >= Max_value){
					Max_value = V[i][j]; Max_r = i; Max_s = j;
				}
			}
		}
	}
	
	// Print the maximum score 
	
	if(type) cout << "\n\nThe maximum alignment score is " << V[0][0] << "\n\n"; // for global alignment 
	else cout << "\n\nThe maximum alignment score is " << Max_value << "\n\n";// for local alignment 
	
	// Output the alignment result
	
	if(type){ // for global alignment
		int i = 0; int j = 0; int a = 0;
		while(i!=L1 && j!=L2){
			if(V[i][j] == V[i+1][j+1]+ (seq1[i]==seq2[j])*match + (seq1[i]!=seq2[j])*mismatch){ 
				align1[a] = seq1[i]; align2[a] = seq2[j]; a++; i++;j++;
			} else{
				if(V[i][j] == V[i+1][j] + indel){ 
					align1[a] = seq1[i]; align2[a] = '_'; a++; i++;
				}
				else {
					align1[a] = '_'; align2[a] = seq2[j]; a++; j++;
				}
			}
		}
		if(i==L1){
			while(j!=L2){
				align1[a] = '_'; align2[a] = seq2[j]; a++;j++;
			}
		} else{
			if(j == L2){
				while(i!=L1){
					align1[a] = seq1[i]; align2[a] = '_'; a++;i++;
				}
			}
		}
		// output the global alignment result
		for(int t = 0; t<=a-1;t++) cout << align1[t];
		cout << "\n";
		for(int t = 0; t<=a-1;t++) cout << align2[t];
		cout << "\n";
	}else{// for local alignment
		int i = Max_r; int j = Max_s; int a = 0; bool cont = true;
		while(cont){
			if(i==L1 || j == L2) break; 
			if((seq1[i]==seq2[j] && V[i][j] == V[i+1][j+1]+ match)||(seq1[i]!=seq2[j] && V[i][j] == V[i+1][j+1]+ mismatch)){
				align1[a] = seq1[i]; align2[a] = seq2[j]; a++; i++;j++;
			} else{
				if(V[i][j] == V[i+1][j] + indel){
					align1[a] = seq1[i]; align2[a] = '_'; a++; i++;
				} else{
					if(V[i][j] == V[i][j+1] + indel){
						align1[a] = '_'; align2[a] = seq2[j]; a++; j++;
					}
					else cont = false; 
				}
			}
		}
		// output the alignment result
		if((Max_r<10 && Max_s >=10 && Max_s <100)||(10<=Max_r && Max_r<100 && Max_s >100)) cout << " "; // for better display
		if((Max_r<10 && Max_s >100)) cout << "  "; // for better display
		cout << Max_r+1 << " ";
		for(int t = 0; t<=a-1;t++) cout << align1[t];
		cout << " " << Max_r + a << "\n";
		if((Max_s<10 && Max_r >=10 && Max_r <100)||(10<=Max_s && Max_s<100 && Max_r >100)) cout << " "; // for better display
		if((Max_s<10 && Max_r >100)) cout << "  "; // for better display
		cout << Max_s+1 << " ";
		for(int t = 0; t<=a-1;t++) cout << align2[t];
		cout << " " << Max_s + a << "\n";
	}
	delete []align1; delete []align2;
	
	cout << "\nPress enter to exit...";
	getchar();
   	getchar();
	return 0;
}