/*
 *                  BioJava development code
 *
 * This code may be freely distributed and modified under the
 * terms of the GNU Lesser General Public Licence.  This should
 * be distributed with the code.  If you do not have a copy,
 * see:
 *
 *      http://www.gnu.org/copyleft/lesser.html
 *
 * Copyright for this code is held jointly by the individual
 * authors.  These should be listed in @author doc comments.
 *
 * For more information on the BioJava project and its aims,
 * or to join the biojava-l mailing list, visit the home page
 * at:
 *
 *      http://www.biojava.org/
 *
 * Created on May 21, 2006
 *
 */
package org.biojava.nbio.structure.align;

import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.ce.GuiWrapper;
import org.biojava.nbio.structure.align.helper.AlignTools;
import org.biojava.nbio.structure.align.helper.JointFragments;
import org.biojava.nbio.structure.align.pairwise.*;
import org.biojava.nbio.structure.geometry.Matrices;
import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.io.PDBFileParser;
import org.biojava.nbio.structure.io.PDBFileReader;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.io.FileOutputStream;
import java.io.InputStream;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.List;

import javax.vecmath.Matrix4d;

/**
 * Perform a pairwise protein structure superimposition.
 *
 * <p>
 * The algorithm is a distance matrix based, rigid body protein structure
 * superimposition. It is based on a variation of the PSC++ algorithm provided
 * by Peter Lackner (Peter.Lackner@sbg.ac.at, personal communication) .
 * </p>
 *
 *
 *
 * <h2>Example</h2>
 * 
 * <pre>
 *  public void run(){
 * 
 * 		// first load two example structures
 * 		{@link InputStream} inStream1 = this.getClass().getResourceAsStream("/files/5pti.pdb");
 * 		{@link InputStream} inStream2 = this.getClass().getResourceAsStream("/files/1tap.pdb");
 * 
 * 		{@link Structure} structure1 = null;
 * 		{@link Structure} structure2 = null;
 * 
 * 		{@link PDBFileParser} pdbpars = new {@link PDBFileParser}();
 * 		structure1 = pdbpars.parsePDBFile(inStream1) ;
 * 		structure2 = pdbpars.parsePDBFile(inStream2);
 * 
 * 
 * 		// calculate structure superimposition for two complete structures
 * 		{@link StructurePairAligner} aligner = new {@link StructurePairAligner}();
 * 
 * 
 * 			// align the full 2 structures with default parameters.
 * 			// see StructurePairAligner for more options and how to align
 * 			// any set of Atoms
 * 			aligner.align(structure1,structure2);
 * 
 * 			{@link AlternativeAlignment}[] aligs = aligner.getAlignments();
 * 			{@link AlternativeAlignment} a = aligs[0];
 * 			System.out.println(a);
 * 
 * 			//display the alignment in Jmol
 * 
 * 			// first get an artificial structure for the alignment
 * 			{@link Structure} artificial = a.getAlignedStructure(structure1, structure2);
 * 
 * 
 * 			// and then send it to Jmol (only will work if Jmol is in the Classpath)
 * 
 * 			BiojavaJmol jmol = new BiojavaJmol();
 * 			jmol.setTitle(artificial.getName());
 * 			jmol.setStructure(artificial);
 * 
 * 			// color the two structures
 * 
 * 
 * 			jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " +
 * 					"select not protein and not solvent; spacefill on;");
 * 			jmol.evalString("select *"+"/1 ; color red; model 1; ");
 * 
 * 
 * 			// now color the equivalent residues ...
 * 
 * 			String[] pdb1 = a.getPDBresnum1();
 * 			for (String res : pdb1 ){
 * 				jmol.evalString("select " + res + "/1 ; backbone 0.6; color white;");
 * 			}
 * 
 * 			jmol.evalString("select *"+"/2; color blue; model 2;");
 * 			String[] pdb2 = a.getPDBresnum2();
 * 			for (String res :pdb2 ){
 * 				jmol.evalString("select " + res + "/2 ; backbone 0.6; color yellow;");
 * 			}
 * 
 * 
 * 			// now show both models again.
 * 			jmol.evalString("model 0;");
 * 
 * 	}
 * </pre>
 *
 *
 *
 * @author Andreas Prlic
 * @author Peter Lackner
 * @since 1.4
 * @version %I% %G%
 */
public class StructurePairAligner {

	private final static Logger logger = LoggerFactory
			.getLogger(StructurePairAligner.class);

	AlternativeAlignment[] alts;
	Matrix distanceMatrix;
	StrucAligParameters params;
	FragmentPair[] fragPairs;

	List<AlignmentProgressListener> listeners = new ArrayList<AlignmentProgressListener>();

	public StructurePairAligner() {
		super();
		params = StrucAligParameters.getDefaultParameters();
		reset();
		alts = new AlternativeAlignment[0];
		distanceMatrix = new Matrix(0, 0);
	}

	public void addProgressListener(AlignmentProgressListener li) {
		listeners.add(li);
	}

	public void clearListeners() {
		listeners.clear();
	}

	/**
	 * example usage of this class
	 *
	 * @param args
	 */
	public static void main(String[] args) throws Exception {
		// UPDATE THE FOLLOWING LINES TO MATCH YOUR SETUP

		PDBFileReader pdbr = new PDBFileReader();
		pdbr.setPath("/Users/andreas/WORK/PDB/");

		// String pdb1 = "1crl";
		// String pdb2 = "1ede";

		String pdb1 = "1buz";
		String pdb2 = "1ali";
		String outputfile = "/tmp/alig_" + pdb1 + "_" + pdb2 + ".pdb";

		// NO NEED TO DO CHANGE ANYTHING BELOW HERE...

		StructurePairAligner sc = new StructurePairAligner();

		// step1 : read molecules

		logger.info("aligning {} vs. {}", pdb1, pdb2);

		Structure s1 = pdbr.getStructureById(pdb1);
		Structure s2 = pdbr.getStructureById(pdb2);

		// step 2 : do the calculations
		sc.align(s1, s2);

		AlternativeAlignment[] aligs = sc.getAlignments();

		// cluster similar results together
		ClusterAltAligs.cluster(aligs);

		// print the result:
		// the AlternativeAlignment object gives also access to rotation
		// matrices / shift vectors.
		for (AlternativeAlignment aa : aligs) {
			logger.info("Alternative Alignment: ", aa);
		}

		// convert AlternativeAlignemnt 1 to PDB file, so it can be opened with
		// a viewer (e.g. Jmol, Rasmol)

		if (aligs.length > 0) {
			AlternativeAlignment aa1 = aligs[0];
			String pdbstr = aa1.toPDB(s1, s2);

			logger.info("writing alignment to {}", outputfile);
			FileOutputStream out = new FileOutputStream(outputfile);
			PrintStream p = new PrintStream(out);

			p.println(pdbstr);

			p.close();
			out.close();
		}

		// display the alignment in Jmol
		// only will work if Jmol is in the Classpath

		if (aligs.length > 0) {

			if (!GuiWrapper.isGuiModuleInstalled()) {
				logger.error("Could not find structure-gui modules in classpath, please install first!");
			}

		}

	}

	private void reset() {
		alts = new AlternativeAlignment[0];
		distanceMatrix = new Matrix(0, 0);
		fragPairs = new FragmentPair[0];

	}

	/**
	 * get the results of step 1 - the FragmentPairs used for seeding the
	 * alignment
	 * 
	 * @return a FragmentPair[] array
	 */

	public FragmentPair[] getFragmentPairs() {
		return fragPairs;
	}

	public void setFragmentPairs(FragmentPair[] fragPairs) {
		this.fragPairs = fragPairs;
	}

	/**
	 * return the alternative alignments that can be found for the two
	 * structures
	 *
	 * @return AlternativeAlignment[] array
	 */
	public AlternativeAlignment[] getAlignments() {
		return alts;
	}

	/**
	 * return the difference of distance matrix between the two structures
	 *
	 * @return a Matrix
	 */
	public Matrix getDistMat() {
		return distanceMatrix;
	}

	/**
	 * get the parameters.
	 *
	 * @return the Parameters.
	 */
	public StrucAligParameters getParams() {
		return params;
	}

	/**
	 * set the parameters to be used for the algorithm
	 *
	 * @param params
	 *            the Parameter object
	 */
	public void setParams(StrucAligParameters params) {
		this.params = params;
	}

	/**
	 * check if debug mode is set on
	 *
	 * @return debug flag
	 */
	@Deprecated
	public boolean isDebug() {
		return false;
	}

	/**
	 * set the debug flag
	 *
	 * @param debug
	 *            flag
	 */
	@Deprecated
	public void setDebug(boolean debug) {
	}

	/**
	 * Calculate the alignment between the two full structures with default
	 * parameters
	 *
	 * @param s1
	 * @param s2
	 * @throws StructureException
	 */
	public void align(Structure s1, Structure s2) throws StructureException {

		align(s1, s2, params);
	}

	/**
	 * Calculate the alignment between the two full structures with user
	 * provided parameters
	 *
	 * @param s1
	 * @param s2
	 * @param params
	 * @throws StructureException
	 */
	public void align(Structure s1, Structure s2, StrucAligParameters params)
			throws StructureException {
		// step 1 convert the structures to Atom Arrays

		Atom[] ca1 = getAlignmentAtoms(s1);
		Atom[] ca2 = getAlignmentAtoms(s2);

		notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
		align(ca1, ca2, params);
	}

	/**
	 * Align two chains from the structures. Uses default parameters.
	 *
	 * @param s1
	 * @param chainId1
	 * @param s2
	 * @param chainId2
	 */
	public void align(Structure s1, String chainId1, Structure s2,
			String chainId2) throws StructureException {
		align(s1, chainId1, s2, chainId2, params);
	}

	/**
	 * Aligns two chains from the structures using user provided parameters.
	 *
	 * @param s1
	 * @param chainId1
	 * @param s2
	 * @param chainId2
	 * @param params
	 * @throws StructureException
	 */
	public void align(Structure s1, String chainId1, Structure s2,
			String chainId2, StrucAligParameters params)
			throws StructureException {
		reset();
		this.params = params;

		Chain c1 = s1.getPolyChainByPDB(chainId1);
		Chain c2 = s2.getPolyChainByPDB(chainId2);

		Structure s3 = new StructureImpl();
		s3.addChain(c1);

		Structure s4 = new StructureImpl();
		s4.addChain(c2);

		Atom[] ca1 = getAlignmentAtoms(s3);
		Atom[] ca2 = getAlignmentAtoms(s4);

		notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
		align(ca1, ca2, params);
	}

	/**
	 * Returns the atoms that are being used for the alignment. (E.g. Calpha
	 * only, etc.)
	 *
	 * @param s
	 * @return an array of Atoms objects
	 */
	public Atom[] getAlignmentAtoms(Structure s) {
		String[] atomNames = params.getUsedAtomNames();
		return StructureTools.getAtomArray(s, atomNames);
	}

	/**
	 * calculate the protein structure superimposition, between two sets of
	 * atoms.
	 *
	 *
	 *
	 * @param ca1
	 *            set of Atoms of structure 1
	 * @param ca2
	 *            set of Atoms of structure 2
	 * @param params
	 *            the parameters to use for the alignment
	 * @throws StructureException
	 */
	public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
			throws StructureException {

		reset();
		this.params = params;

		long timeStart = System.currentTimeMillis();

		// step 1 get all Diagonals of length X that are similar between both
		// structures
		logger.debug(" length atoms1:" + ca1.length);
		logger.debug(" length atoms2:" + ca2.length);

		logger.debug("step 1 - get fragments with similar intramolecular distances ");

		int k = params.getDiagonalDistance();
		int k2 = params.getDiagonalDistance2();
		int fragmentLength = params.getFragmentLength();

		if (ca1.length < (fragmentLength + 1)) {
			throw new StructureException("structure 1 too short (" + ca1.length
					+ "), can not align");
		}
		if (ca2.length < (fragmentLength + 1)) {
			throw new StructureException("structure 2 too short (" + ca2.length
					+ "), can not align");
		}
		int rows = ca1.length - fragmentLength + 1;
		int cols = ca2.length - fragmentLength + 1;
		distanceMatrix = new Matrix(rows, cols, 0.0);

		double[] dist1 = AlignTools.getDiagonalAtK(ca1, k);

		double[] dist2 = AlignTools.getDiagonalAtK(ca2, k);
		double[] dist3 = new double[0];
		double[] dist4 = new double[0];
		if (k2 > 0) {
			dist3 = AlignTools.getDiagonalAtK(ca1, k2);
			dist4 = AlignTools.getDiagonalAtK(ca2, k2);
		}

		double[][] utmp = new double[][] { { 0, 0, 1 } };
		Atom unitvector = new AtomImpl();
		unitvector.setCoords(utmp[0]);

		List<FragmentPair> fragments = new ArrayList<FragmentPair>();

		for (int i = 0; i < rows; i++) {

			Atom[] catmp1 = AlignTools.getFragment(ca1, i, fragmentLength);
			Atom center1 = AlignTools.getCenter(ca1, i, fragmentLength);

			for (int j = 0; j < cols; j++) {

				double rdd1 = AlignTools.rms_dk_diag(dist1, dist2, i, j,
						fragmentLength, k);
				double rdd2 = 0;
				if (k2 > 0)
					rdd2 = AlignTools.rms_dk_diag(dist3, dist4, i, j,
							fragmentLength, k2);
				double rdd = rdd1 + rdd2;
				distanceMatrix.set(i, j, rdd);

				if (rdd < params.getFragmentMiniDistance()) {
					FragmentPair f = new FragmentPair(fragmentLength, i, j);
					Atom[] catmp2 = AlignTools.getFragment(ca2, j,
							fragmentLength);
					Atom center2 = AlignTools.getCenter(ca2, j,
							fragmentLength);

					f.setCenter1(center1);
					f.setCenter2(center2);

					Matrix4d t = SuperPositions.superpose(
							Calc.atomsToPoints(catmp1),
							Calc.atomsToPoints(catmp2));
					
					Matrix rotmat = Matrices.getRotationJAMA(t);
					f.setRot(rotmat);

					Atom aunitv = (Atom) unitvector.clone();
					Calc.rotate(aunitv, rotmat);
					f.setUnitv(aunitv);

					boolean doNotAdd = false;
					if (params.reduceInitialFragments()) {
						doNotAdd = FragmentJoiner.reduceFragments(
								fragments, f, distanceMatrix);

					}
					if (doNotAdd)
						continue;

					fragments.add(f);
				}
			}
		}

		notifyFragmentListeners(fragments);

		FragmentPair[] fp = fragments
				.toArray(new FragmentPair[fragments.size()]);
		setFragmentPairs(fp);

		logger.debug(" got # fragment pairs: {}", fp.length);

		logger.debug("step 2 - join fragments");

		// step 2 combine them to possible models
		FragmentJoiner joiner = new FragmentJoiner();

		JointFragments[] frags;

		if (params.isJoinFast()) {
			// apply the quick alignment procedure.
			// less quality in alignments, better for DB searches...
			frags = joiner.approach_ap3(ca1, ca2, fp, params);

			joiner.extendFragments(ca1, ca2, frags, params);

		} else if (params.isJoinPlo()) {
			// this approach by StrComPy (peter lackner):
			frags = joiner.frag_pairwise_compat(fp, params.getAngleDiff(),
					params.getFragCompat(), params.getMaxrefine());

		} else {

			// my first implementation
			frags = joiner.approach_ap3(ca1, ca2, fp, params);
		}

		notifyJointFragments(frags);

		logger.debug(" number joint fragments: ", frags.length);

		logger.debug("step 3 - refine alignments");

		List<AlternativeAlignment> aas = new ArrayList<AlternativeAlignment>();
		for (int i = 0; i < frags.length; i++) {
			JointFragments f = frags[i];
			AlternativeAlignment a = new AlternativeAlignment();
			a.apairs_from_idxlst(f);
			a.setAltAligNumber(i + 1);
			a.setDistanceMatrix(distanceMatrix);

			try {
				if (params.getMaxIter() > 0) {

					a.refine(params, ca1, ca2);
				} else {

					a.finish(params, ca1, ca2);

				}
			} catch (StructureException e) {
				logger.error("Refinement of fragment {} failed", i, e);
			}
			a.calcScores(ca1, ca2);
			aas.add(a);
		}

		// sort the alternative alignments
		Comparator<AlternativeAlignment> comp = new AltAligComparator();
		Collections.sort(aas, comp);
		Collections.reverse(aas);

		alts = aas.toArray(new AlternativeAlignment[aas.size()]);
		// do final numbering of alternative solutions
		int aanbr = 0;
		for (AlternativeAlignment a : alts) {
			aanbr++;
			a.setAltAligNumber(aanbr);
		}

		logger.debug("total calculation time: {} ms.",
				(System.currentTimeMillis() - timeStart));
	}

	private void notifyStartingAlignment(String name1, Atom[] ca1,
			String name2, Atom[] ca2) {
		for (AlignmentProgressListener li : listeners) {
			li.startingAlignment(name1, ca1, name2, ca2);
		}
	}

	private void notifyFragmentListeners(List<FragmentPair> fragments) {

		for (AlignmentProgressListener li : listeners) {
			li.calculatedFragmentPairs(fragments);
		}

	}

	private void notifyJointFragments(JointFragments[] fragments) {
		for (AlignmentProgressListener li : listeners) {
			li.jointFragments(fragments);
		}
	}

}
