/*
 * This file is part of Jstacs.
 * 
 * Jstacs is free software: you can redistribute it and/or modify it under the
 * terms of the GNU General Public License as published by the Free Software
 * Foundation, either version 3 of the License, or (at your option) any later
 * version.
 * 
 * Jstacs is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 * A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License along with
 * Jstacs. If not, see <http://www.gnu.org/licenses/>.
 * 
 * For more information on Jstacs, visit http://www.jstacs.de
 */

package de.jstacs.algorithms;

import de.jstacs.algorithms.Alignment.Costs.Direction;
import de.jstacs.data.AlphabetContainer;
import de.jstacs.data.Sequence;

/**
 * Class for gapped global alignments using the extended Needleman-Wunsch
 * algorithm.
 * 
 * @author Jan Grau
 * 
 */
public class Alignment {

	private Sequence s1, s2;

	private Costs costs;

	private D[][] d;

	private E[][] e;

	private F[][] f;

	/**
	 * Creates a new {@link Alignment} instance that aligns the sequences
	 * <code>s1</code> and <code>s2</code> using the costs defined in
	 * <code>costs</code>.
	 * 
	 * @param s1
	 *            the first sequence
	 * @param s2
	 *            the second sequence
	 * @param costs
	 *            the costs
	 */
	public Alignment( Sequence s1, Sequence s2, Costs costs ) {
		this.s1 = s1;
		this.s2 = s2;
		this.costs = costs;
	}

	/**
	 * Computes and returns the alignment of <code>s1</code> and <code>s2</code>
	 * ({@link #Alignment(Sequence, Sequence, Costs)}).
	 * 
	 * @return the alignment
	 */
	public StringAlignment getAlignment() {

		d = new D[s1.getLength() + 1][s2.getLength() + 1];
		e = new E[s1.getLength() + 1][s2.getLength() + 1];
		f = new F[s1.getLength() + 1][s2.getLength() + 1];

		for( int i = 0; i <= s1.getLength(); i++ ) {
			for( int j = 0; j <= s2.getLength(); j++ ) {
				e[i][j] = new E( i, j );
				f[i][j] = new F( i, j );
				d[i][j] = new D( i, j );
			}
		}

		Element curr = null;
		if( e[s1.getLength()][s2.getLength()].cost < f[s1.getLength()][s2.getLength()].cost && e[s1.getLength()][s2.getLength()].cost < d[s1.getLength()][s2.getLength()].cost ) {
			curr = e[s1.getLength()][s2.getLength()];
		} else if( f[s1.getLength()][s2.getLength()].cost < d[s1.getLength()][s2.getLength()].cost ) {
			curr = f[s1.getLength()][s2.getLength()];
		} else {
			curr = d[s1.getLength()][s2.getLength()];
		}

		double cost = curr.cost;

		StringBuffer b1 = new StringBuffer();
		StringBuffer b2 = new StringBuffer();

		AlphabetContainer cont = s1.getAlphabetContainer();

		while( curr.pre != null ) {
			if( curr.direction == Direction.DIAGONAL ) {
				b1.insert( 0, cont.getSymbol( curr.i - 1, s1.discreteVal( curr.i - 1 ) ) );
				b2.insert( 0, cont.getSymbol( curr.j - 1, s2.discreteVal( curr.j - 1 ) ) );
			} else if( curr.direction == Direction.LEFT ) {
				b1.insert( 0, '-' );
				b2.insert( 0, cont.getSymbol( curr.j - 1, s2.discreteVal( curr.j - 1 ) ) );
			} else if( curr.direction == Direction.TOP ) {
				b1.insert( 0, cont.getSymbol( curr.j - 1, s2.discreteVal( curr.j - 1 ) ) );
				b2.insert( 0, '-' );
			}
			curr = curr.pre;
		}

		return new StringAlignment( b1.toString(), b2.toString(), cost );

	}

	/**
	 * Class for the representation of an alignment of two {@link String}s. It
	 * contains the two {@link String}s that were aligned and expanded by
	 * gap-symbols, and the edit-costs according to the employed
	 * {@link Alignment.Costs} instance.
	 * 
	 * @author Jan Grau
	 */
	public static class StringAlignment {

		private String r1;

		private String r2;

		private double cost;

		/**
		 * Creates the instance for the two (extended) {@link String}s and the
		 * edit-costs.
		 * 
		 * @param r1
		 *            the first {@link String}
		 * @param r2
		 *            the second {@link String}
		 * @param cost
		 *            the edit-costs
		 */
		protected StringAlignment( String r1, String r2, double cost ) {
			this.r1 = r1;
			this.r2 = r2;
			this.cost = cost;
		}

		/**
		 * Returns the costs.
		 * 
		 * @return the costs
		 */
		public double getCost() {
			return cost;
		}

		/**
		 * Returns the first {@link String}.
		 * 
		 * @return the first {@link String}
		 */
		public String getFirst() {
			return r1;
		}

		/**
		 * Returns the second {@link String}.
		 * 
		 * @return the second {@link String}
		 */
		public String getSecond() {
			return r2;
		}

		/* (non-Javadoc)
		 * @see java.lang.Object#toString()
		 */
		@Override
		public String toString() {
			StringBuffer buf = new StringBuffer( r1 );
			buf.append( "\n" );
			buf.append( r2 );
			buf.append( "\n" );
			buf.append( "costs: " );
			buf.append( cost );
			return buf.toString();
		}

	}

	/**
	 * Class for a general element of the DP (<b>d</b>ynamic
	 * <b>p</b>rogramming)-matrices of the Needleman-Wunsch algorithm.
	 * 
	 * @author Jan Grau
	 */
	private abstract class Element {

		/**
		 * The row index of the element in the matrix.
		 */
		protected int i;

		/**
		 * The column index of the element in the matrix.
		 */
		protected int j;

		/**
		 * The costs of the alignment until this element, i.e.
		 * <code>s1(1..i),s2(1..j)</code>.
		 */
		double cost;

		/**
		 * The direction of the predecessor in the DP-matrix.
		 */
		protected Direction direction;

		/**
		 * The predecessor-element in the DP-matrix.
		 */
		protected Element pre;

		/**
		 * Creates a new element for row <code>i</code> and column
		 * <code>j</code>.
		 * 
		 * @param i
		 *            the row index for the new element
		 * @param j
		 *            the column index for the new element
		 */
		public Element( int i, int j ) {
			this.i = i;
			this.j = j;
		}

	}

	private class D extends Element {

		private D( int i, int j ) {
			super( i, j );
			if( i == 0 && j == 0 ) {
				cost = 0;
			} else if( i == 0 && j > 0 ) {
				cost = Double.POSITIVE_INFINITY;
				direction = Direction.LEFT;
			} else if( i > 0 && j == 0 ) {
				cost = Double.POSITIVE_INFINITY;
				direction = Direction.TOP;
			} else {
				double diag = d[i - 1][j - 1].cost + costs.getCostFor( s1, s2, i, j, Direction.DIAGONAL );
				double left = e[i][j].cost;
				double top = f[i][j].cost;

				if( diag < left && diag < top ) {
					cost = diag;
					direction = Direction.DIAGONAL;
					pre = d[i - 1][j - 1];
				} else if( left < top ) {
					cost = left;
					direction = Direction.SELF;
					pre = e[i][j];
				} else {
					cost = top;
					direction = Direction.SELF;
					pre = f[i][j];
				}
			}

		}

	}

	private class E extends Element {

		private E( int i, int j ) {
			super( i, j );
			if( j == 0 ) {
				cost = Double.POSITIVE_INFINITY;
			} else if( i == 0 && j > 0 ) {
				direction = Direction.LEFT;
				cost = costs.getGapCostsFor( j );
				pre = e[i][j - 1];
			} else {
				double elong = e[i][j - 1].cost + costs.getElongateCosts();
				double start = d[i][j - 1].cost + costs.getGapCostsFor( 1 );
				if( elong < start ) {
					cost = elong;
					direction = Direction.LEFT;
					pre = e[i][j - 1];
				} else {
					cost = start;
					direction = Direction.LEFT;
					pre = d[i][j - 1];
				}
			}
		}

	}

	private class F extends Element {

		private F( int i, int j ) {
			super( i, j );
			if( i == 0 ) {
				cost = Double.POSITIVE_INFINITY;
			} else if( i > 0 && j == 0 ) {
				direction = Direction.TOP;
				cost = costs.getGapCostsFor( j );
				pre = f[i - 1][j];
			} else {
				double elong = f[i - 1][j].cost + costs.getElongateCosts();
				double start = d[i - 1][j].cost + costs.getGapCostsFor( 1 );
				if( elong < start ) {
					cost = elong;
					direction = Direction.TOP;
					pre = f[i - 1][j];
				} else {
					cost = start;
					direction = Direction.TOP;
					pre = d[i - 1][j];
				}
			}
		}

	}

	/**
	 * Class for simple costs with costs <code>mismatch</code> for a mismatch,
	 * costs <code>start</code> to start a new gap, costs <code>elong</code> to
	 * elongate a gap by one position and costs of <code>0</code> for a match.
	 */
	public static class SimpleCosts implements Costs {

		private double mismatch;

		private double start;

		private double elong;

		/**
		 * Creates a new instance of simple costs with costs
		 * <code>mismatch</code> for a mismatch, costs <code>start</code> to
		 * start a new gap, costs <code>elong</code> to elongate a gap by one
		 * position and costs of <code>0</code> for a match.
		 * 
		 * @param mismatch
		 *            the mismatch costs
		 * @param start
		 *            the costs to start a gap
		 * @param elong
		 *            the costs to elongate a gap
		 */
		public SimpleCosts( double mismatch, double start, double elong ) {
			this.start = start;
			this.elong = elong;
			this.mismatch = mismatch;
		}

		/* (non-Javadoc)
		 * @see de.jstacs.algorithms.Alignment.Costs#getCostFor(de.jstacs.data.Sequence, de.jstacs.data.Sequence, int, int, de.jstacs.algorithms.Alignment.Costs.Direction)
		 */
		public double getCostFor( Sequence s1, Sequence s2, int i, int j, Direction from ) {
			if( from == Direction.TOP || from == Direction.LEFT ) {
				return start + elong;
			} else {
				if( s1.discreteVal( i - 1 ) != s2.discreteVal( j - 1 ) ) {
					return mismatch;
				} else {
					return 0;
				}
			}
		}

		/* (non-Javadoc)
		 * @see de.jstacs.algorithms.Alignment.Costs#getElongateCosts()
		 */
		public double getElongateCosts() {
			return elong;
		}

		/* (non-Javadoc)
		 * @see de.jstacs.algorithms.Alignment.Costs#getGapCostsFor(int)
		 */
		public double getGapCostsFor( int length ) {
			return start + ( length * elong );
		}

	}

	/**
	 * The general interface for the costs of an alignment.
	 * 
	 * @author Jan Grau
	 */
	public static interface Costs {

		/**
		 * The direction of the predecessor in the DP (<b>d</b>ynamic
		 * <b>p</b>rogramming)-matrix.
		 * 
		 * @author Jan Grau
		 */
		public enum Direction {

			/**
			 * The predecessor is above the current element.
			 */
			TOP,
			/**
			 * The predecessor is left of the current element.
			 */
			LEFT,
			/**
			 * The predecessor is top left of the current element.
			 */
			DIAGONAL,
			/**
			 * The predecessor is at the same position as the current element.
			 */
			SELF;

		}

		/**
		 * Returns the costs for the alignment if <code>s1(i)</code> and
		 * <code>s2(j)</code> coming from <code>from</code>.
		 * 
		 * @param s1
		 *            the first sequence
		 * @param s2
		 *            the second sequence
		 * @param i
		 *            the index in the first sequence
		 * @param j
		 *            the index in the second sequence
		 * @param from
		 *            the direction within the DP-matrix
		 * 
		 * @return the costs
		 */
		public double getCostFor( Sequence s1, Sequence s2, int i, int j, Direction from );

		/**
		 * Returns the costs for a gap of length <code>length</code>.
		 * 
		 * @param length
		 *            the length of the gap
		 * 
		 * @return the corresponding costs
		 */
		public double getGapCostsFor( int length );

		/**
		 * Returns the costs to elongate a gap by one position.
		 * 
		 * @return the corresponding costs
		 */
		public double getElongateCosts();

	}

}
