/*******************************************************************************
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * This library 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 Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 * 
 * Contact Info:
 * 	Bruce Donald
 * 	Duke University
 * 	Department of Computer Science
 * 	Levine Science Research Center (LSRC)
 * 	Durham
 * 	NC 27708-0129 
 * 	USA
 * 	brd@cs.duke.edu
 * 
 * Copyright (C) 2011 Jeffrey W. Martin and Bruce R. Donald
 * 
 * <signature of Bruce Donald>, April 2011
 * Bruce Donald, Professor of Computer Science
 ******************************************************************************/


package edu.duke.donaldLab.jdshot.disco;

import java.util.ArrayList;
import java.util.Collections;
import java.util.PriorityQueue;

import Jama.EigenvalueDecomposition;
import Jama.Matrix;
import edu.duke.donaldLab.share.geom.Vector3;
import edu.duke.donaldLab.share.kinemage.Kinemage;
import edu.duke.donaldLab.share.kinemage.KinemageBuilder;
import edu.duke.donaldLab.share.kinemage.KinemageWriter;
import edu.duke.donaldLab.share.math.EigPair;
import edu.duke.donaldLab.share.protein.AtomAddressInternal;
import edu.duke.donaldLab.share.protein.Protein;
import edu.duke.donaldLab.share.protein.Subunit;
import edu.duke.donaldLab.share.protein.tools.ProteinGeometry;

public class ReferenceAnalyzer
{
	/**************************
	 *   Static Methods
	 **************************/
	
	public static Vector3 getSymmetryAxisPosition( Protein oligomer )
	{
		return ProteinGeometry.getCentroid( oligomer, oligomer.getBackboneAtomIndex() );
	}
	
	public static Vector3 getSymmetryAxisOrientation( Protein oligomer )
	{
		// work with a copy of the protein so we can modify it
		oligomer = new Protein( oligomer );
		
		// center the protein on the average backbone position (aka zero mean)
		ProteinGeometry.center( oligomer, oligomer.getBackboneAtomIndex() );
		
		// average the positions of all symmetric backbone atoms
		ArrayList<Vector3> atoms = new ArrayList<Vector3>( oligomer.getSubunit( 0 ).getBackboneAtomIndex().size() );
		for( AtomAddressInternal address : oligomer.getSubunit( 0 ).getBackboneAtomIndex() )
		{
			Vector3 tempPos = new Vector3( oligomer.getAtom( address ).getPosition() );
			for( int i=1; i<oligomer.getSubunits().size(); i++ )
			{
				tempPos.add( oligomer.getAtom(
					i,
					address.getResidueId(),
					address.getAtomId()
				).getPosition() );
			}
			tempPos.scale( 1.0 / oligomer.getSubunits().size() );
			atoms.add( tempPos );
		}
		
		// do PCA on the backbone atom positions
		Matrix m = new Matrix( 3, atoms.size() );
		for( int i=0; i<atoms.size(); i++ )
		{
			Vector3 pos = atoms.get( i );
			m.set( 0, i, pos.x );
			m.set( 1, i, pos.y );
			m.set( 2, i, pos.z );
		}
		Matrix cor = m.times( m.transpose() );
		EigenvalueDecomposition eig = new EigenvalueDecomposition( cor );
		
		// find the largest and second-largest eigenvalues/eigenvectors
		PriorityQueue<EigPair> q = new PriorityQueue<EigPair>( 3, Collections.reverseOrder() );
		q.add( new EigPair( eig, 0 ) );
		q.add( new EigPair( eig, 1 ) );
		q.add( new EigPair( eig, 2 ) );
		EigPair z = q.poll();
		EigPair y = q.poll();
		
		// TEMP: show a kinemage of the axis, protein, and positions
		Kinemage kinemage = new Kinemage();
		KinemageBuilder.appendBackbone( kinemage, oligomer );
		KinemageBuilder.appendPoints( kinemage, atoms, "Points", 3, 6 );
		KinemageBuilder.appendVector(
			kinemage,
			z.getEigenvector(),
			ProteinGeometry.getCentroid( oligomer ),
			"Axis",
			1,
			3,
			10
		);
		new KinemageWriter().showAndWait( kinemage );
		
		// TEMP
		System.out.println( z.getEigenvalue() / y.getEigenvalue() );
		
		// check the differences between eigenvalues
		// only return confident orientations
		if( z.getEigenvalue() / y.getEigenvalue() < 10.0 )
		{
			return null;
		}
		
		return z.getEigenvector();
	}
	
	public static Vector3 getDnSymmetryAxisOrientation( Protein oligomer, Vector3 cnAxis )
	{
		// look at all pairs of subunits that include subunit 0
		for( int i=1; i<oligomer.getSubunits().size(); i++ )
		{
			Subunit a = oligomer.getSubunit( 0 );
			Subunit b = oligomer.getSubunit( i );
			
			// build a temp protein with subunits a and b
			ArrayList<Subunit> subunits  = new ArrayList<Subunit>( 2 );
			subunits.add( a );
			subunits.add( b );
			Protein protein = new Protein();
			protein.setSubunits( subunits );
			
			// UNDONE: not sure what this will return...
			// THIS DOESN'T WORK.
			Vector3 axis = getSymmetryAxisOrientation( protein );
			
			// is it parallel to the cn axis?
			if( Math.abs( axis.getDot( cnAxis ) ) > 0.9 )
			{
				continue;
			}
			
			return axis;
		}
		
		return null;
	}
	
	public Kinemage visualizeSymmetryAxis( Protein protein )
	throws Exception
	{
		Kinemage kinemage = new Kinemage( "Symmetry Axis Visualization" );
		KinemageBuilder.appendAxes( kinemage );
		KinemageBuilder.appendBackbone(
			kinemage,
			protein,
			protein.getName(),
			0,
			2
		);
		KinemageBuilder.appendVector(
			kinemage,
			getSymmetryAxisPosition( protein ),
			getSymmetryAxisOrientation( protein ),
			"Symmetry Axis",
			1,
			3,
			10
		);
		return kinemage;
	}
}
