/*******************************************************************************
 * 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.HashMap;
import java.util.List;
import java.util.Map;

import org.apache.log4j.Logger;

import edu.duke.donaldLab.jdshot.disco.AnnulusIndex.AnnulusEntry;
import edu.duke.donaldLab.share.geom.Annulus2;
import edu.duke.donaldLab.share.geom.Vector2;
import edu.duke.donaldLab.share.geom.Vector3;
import edu.duke.donaldLab.share.io.Logging;
import edu.duke.donaldLab.share.math.Matrix2;
import edu.duke.donaldLab.share.nmr.Assignment;
import edu.duke.donaldLab.share.nmr.DistanceRestraint;
import edu.duke.donaldLab.share.nmr.DistanceRestraintMapper;
import edu.duke.donaldLab.share.nmr.DistanceRestraintWriter;
import edu.duke.donaldLab.share.protein.AtomAddressInternal;
import edu.duke.donaldLab.share.protein.AtomAddressReadable;
import edu.duke.donaldLab.share.protein.Protein;
import edu.duke.donaldLab.share.protein.Subunit;
import edu.duke.donaldLab.share.protein.tools.MonomerCloner;

public class AnnulusCalculator
{
	/**************************
	 *   Definitions
	 **************************/
	
	private static final double Factor = 1e10;
	
	
	/**************************
	 *   Data Members
	 **************************/
	
	private static final Logger m_log = Logging.getLog( AnnulusCalculator.class );
	
	
	/**************************
	 *   Static Methods
	 **************************/
	
	public static AnnulusIndex computeAnnuli( Subunit subunit, List<DistanceRestraint<AtomAddressInternal>> restraints, int numSubunits )
	{
		int numAssignments = 0;
		int numAnnuli = 0;
		AnnulusIndex index = new AnnulusIndex();
		for( DistanceRestraint<AtomAddressInternal> restraint : restraints )
		{
			numAssignments += restraint.getNumAssignments();
			
			// one side of the restraint should reference only one subunit
			int numLeftSubunits = AtomAddressInternal.getNumDistinctSubunits( restraint.getLefts() );
			int numRightSubunits = AtomAddressInternal.getNumDistinctSubunits( restraint.getRights() );
			assert( numLeftSubunits == 1 || numRightSubunits == 1 );
			
			// create a copy of the restraint so we can modify it
			DistanceRestraint<AtomAddressInternal> copy = new DistanceRestraint<AtomAddressInternal>( restraint );
			if( numLeftSubunits != 1 && numRightSubunits == 1 )
			{
				copy.swap();
			}
			assert( AtomAddressInternal.getNumDistinctSubunits( copy.getLefts() ) == 1 );
			
			// rotate the restraint, so the left side only references the left subunit
			int leftSubunitId = copy.getLefts().iterator().next().getSubunitId();
			AtomAddressInternal.rotateSubunits( copy.getLefts(), numSubunits, -leftSubunitId );
			AtomAddressInternal.rotateSubunits( copy.getRights(), numSubunits, -leftSubunitId );
			assert( AtomAddressInternal.getNumDistinctSubunits( copy.getLefts() ) == 1 );
			assert( copy.getLefts().iterator().next().getSubunitId() == 0 );
			
			for( Assignment<AtomAddressInternal> assignment : copy )
			{
				// just ignore intra-molecular assignments
				if( !assignment.isIntermolecular() )
				{
					continue;
				}
				
				// compute the annulus
				Annulus2 annulus = computeAnnulus(
					subunit,
					assignment,
					copy.getMinDistance(),
					copy.getMaxDistance(),
					numSubunits
				);
				index.add( restraint, assignment, annulus );
				numAnnuli++;
			}
		}
		
		// LOGGING
		m_log.info( String.format( "Computed %d annuli from %d assignments", numAnnuli, numAssignments ) );
		
		return index;
	}
	
	public static AnnulusIndex computeUntrustedAnnuli( Subunit subunit, List<DistanceRestraint<AtomAddressInternal>> untrustedRestraints, int numSubunits, List<FaceInfo> msrs )
	{
		AnnulusIndex untrustedAnnulusIndex = computeAnnuli( subunit, untrustedRestraints, numSubunits );
		
		// discard annuli if they don't intersect any of the MSRs
		AnnulusIndex filteredAnnulusIndex = new AnnulusIndex();
		for( Map.Entry<DistanceRestraint<AtomAddressInternal>,HashMap<Assignment<AtomAddressInternal>,AnnulusEntry>> entry : untrustedAnnulusIndex )
		{
			for( Map.Entry<Assignment<AtomAddressInternal>,AnnulusEntry> subEntry : entry.getValue().entrySet() )
			{
				for( FaceInfo msr : msrs )
				{
					if( msr.intersects( subEntry.getValue().annulus ) )
					{
						filteredAnnulusIndex.add( entry.getKey(), subEntry.getKey(), subEntry.getValue().annulus );
						break;
					}
				}
			}
		}
		
		// LOGGING
		m_log.info( String.format( "Untrusted distance restraints consistent with the trusted MSRs: %d/%d (%.1f%%)",
			filteredAnnulusIndex.getNumRestraints(),
			untrustedRestraints.size(),
			(double)filteredAnnulusIndex.getNumRestraints() / (double)untrustedRestraints.size() * 100.0
		) );
		
		return filteredAnnulusIndex;
	}
	
	public static void reportUnsatisfiableRestraints( AnnulusIndex annulusIndex, Subunit subunit, int numSubunits )
	{
		// search for emtpy annuli
		ArrayList<DistanceRestraint<AtomAddressInternal>> badRestraints = new ArrayList<DistanceRestraint<AtomAddressInternal>>();
		for( Map.Entry<DistanceRestraint<AtomAddressInternal>,HashMap<Assignment<AtomAddressInternal>,AnnulusEntry>> entry : annulusIndex )
		{
			boolean isOneAssignmentGood = false;
			for( Map.Entry<Assignment<AtomAddressInternal>,AnnulusEntry> subEntry : entry.getValue().entrySet() )
			{
				if( subEntry.getValue().annulus.maxRadius > 0.0 )
				{
					isOneAssignmentGood = true;
					break;
				}
			}
			
			if( !isOneAssignmentGood )
			{
				badRestraints.add( entry.getKey() );
			}
		}
		
		// report empty annuli
		if( badRestraints.size() > 0 )
		{
			Protein oligomer = MonomerCloner.clone( subunit, numSubunits );
			DistanceRestraintWriter writer = new DistanceRestraintWriter();
			StringBuilder buf = new StringBuilder();
			buf.append( badRestraints.size() );
			buf.append( " distance restraints are unsatisfiable using the computed symmetry axis orientation:\n" );
			for( DistanceRestraint<AtomAddressInternal> restraint : badRestraints )
			{
				DistanceRestraint<AtomAddressReadable> readableRestraint = DistanceRestraintMapper.mapInternalToReadable( restraint, oligomer );
				buf.append( "\t" );
				buf.append( writer.writeToString( readableRestraint ) );
				buf.append( "\n" );
			}
			m_log.warn( buf.toString() );
		}
	}
	

	/**************************
	 *   Static Functions
	 **************************/
	
	private static Annulus2 computeAnnulus( Subunit subunit, Assignment<AtomAddressInternal> assignment, double minDist, double maxDist, int numSubunits )
	{
		// here's the general formula
		// A = ( 1/h W^-1( Rq_A-p' ), r_l/h, r_u/h )
		
		// make sure the first subunit is on the left
		assert( assignment.getLeft().getSubunitId() == 0 );
		int subunitId = assignment.getRight().getSubunitId();
		
		/* get the non-zero subunit
		assert( assignment.getLeft().getSubunitId() == 0 || assignment.getRight().getSubunitId() == 0 );
		int subunitId = assignment.getLeft().getSubunitId();
		if( subunitId == 0 )
		{
			subunitId = assignment.getRight().getSubunitId();
		}
		assert( subunitId != 0 );
		*/
		
		// get the atom positions q_A and p'
		Vector3 p = subunit.getAtom( assignment.getLeft() ).getPosition();
		Vector2 pPrime = new Vector2( p.x, p.y );
		Vector3 qA3 = subunit.getAtom( assignment.getRight() ).getPosition();
		Vector2 qA = new Vector2( qA3.x, qA3.y );
		
		// compute the rotation R
		double alpha = Math.PI * 2.0 / (double)numSubunits;
		Matrix2 R = new Matrix2();
		Matrix2.getRotation( R, alpha * (double)subunitId );
		
		// compute the rotation W
		Matrix2 W = new Matrix2(
			R.data[0][0] - 1.0, R.data[0][1],
			R.data[1][0], R.data[1][1] - 1.0
		);
		
		// compute the scaling h
		double uMag = new Vector2( W.data[0][0], W.data[1][0] ).getLength();
		double vMag = new Vector2( W.data[0][1], W.data[1][1] ).getLength();
		assert( uMag == vMag );
		double h = uMag;
		
		// scale W
		W.data[0][0] /= h;
		W.data[0][1] /= h;
		W.data[1][0] /= h;
		W.data[1][1] /= h;
		
		// compute the center: 1/h W^-1( Rq_A-p' )
		Vector2 center = qA;
		R.multiply( center );
		center.subtract( pPrime );
		W.invert();
		W.multiply( center );
		center.scale( 1.0 / h );
		
		// compute the radii
		double rl = computeRadius( p, qA3.z, minDist ) / h;
		double ru = computeRadius( p, qA3.z, maxDist ) / h;
		
		/* Note:
			slight numerical instability is causing problems with the CGAL arrangement
			package which is sensitive enough to detect this insignificant changes.
			So, round all numbers to 10 decimal places. That should be enough
			precision.
		*/
		center.x = Math.round( center.x * Factor ) / Factor;
		center.y = Math.round( center.y * Factor ) / Factor;
		rl = Math.round( rl * Factor ) / Factor;
		ru = Math.round( ru * Factor ) / Factor;

		return new Annulus2( center, rl, ru );
	}
	
	private static double computeRadius( Vector3 p, double z, double dist )
	{
		z -= p.z;
		if( Math.abs( z ) >= dist )
		{
			return 0.0;
		}
		return Math.sqrt( dist*dist - z*z );
	}	
}
