//  $Id: Octree.java,v 1.4 2007/02/18 05:08:05 Sasha Buzko Exp $
//
//  Copyright (c) 2000-2003  San Diego Supercomputer Center (SDSC),
//  a facility operated jointly by the University of California,
//  San Diego (UCSD) and General Atomics, San Diego, California, USA.
//
//  Users and possessors of this source code are hereby granted a
//  nonexclusive, royalty-free copyright and design patent license to
//  use this code in individual software.  License is not granted for
//  commercial resale, in whole or in part, without prior written
//  permission from SDSC.  This source is provided "AS IS" without express
//  or implied warranty of any kind.
//
//  For further information, please see:  http://mbt.sdsc.edu
//
//  History:
//  $Log: Octree.java,v $
//  Revision 1.4  2007/02/18 05:08:05  Sasha Buzko
//  *** empty log message ***
//
//  Revision 1.3  2007/02/07 18:42:57  Sasha Buzko
//  *** empty log message ***
//
//  Revision 1.2  2007/02/07 02:45:49  Sasha Buzko
//  *** empty log message ***
//
//  Revision 1.1  2006/05/20 17:02:04  Sasha Buzko
//  Updated version
//
//  Revision 1.1  2006/04/30 20:13:59  Sasha Buzko
//  New version of the app
//
//  Revision 1.1  2006/04/15 19:42:26  Sasha Buzko
//  Initial commit
//
//  Revision 1.2  2005/11/26 04:36:10  Administrator
//  *** empty log message ***
//
//  Revision 1.1  2005/11/13 04:35:04  Administrator
//  *** empty log message ***
//
//  Revision 1.5  2003/10/17 18:19:32  moreland
//  Fixed a javadoc comment.
//
//  Revision 1.4  2003/10/01 21:11:51  agramada
//  Added methods needed by DerivedInformation class for derivation of
//  secondary structures.
//
//  Revision 1.3  2003/07/17 16:40:32  agramada
//  Removed some comments.
//
//  Revision 1.2  2003/07/11 20:23:05  agramada
//  Made get*Bonds methods return Bond objects instead of BondInfo. Removed
//  some coments.
//
//  Revision 1.1  2003/07/11 18:17:53  moreland
//  Modifed Apostol's Octree classes to genate Bonds from the BondFactory
//  and in turn the StructureMap class.
//
//  Revision 1.3  2003/07/07 23:02:47  agramada
//  Added a method to return a vector of covalent bonds rather than an array.
//
//  Revision 1.2  2003/06/24 22:19:46  agramada
//  Reorganized the geometry package. Old classes removed, new classes added.
//
//  Revision 1.1  2003/04/24 18:46:36  agramada
//  First version of the Octree class to be used for an efficient search of
//  the set of bonds between a set of atoms.
//
//  Revision 1.2  2003/02/20 17:08:22  agramada
//  First working version.
//
//  Revision 1.0  2003/02/20 23:38:39  agramada
//


package edu.sdsc.mbt.util;


import java.util.*;
import edu.sdsc.mbt.Bond;


/**
 * Octree class constructed with recursion in mind.
 * Each child is itself a tree.
 * <P>
 * @author      Apostol Gramada
 */
public class Octree
{
	private int dimension;  // Keep it general as long as practical. In practice, maximum 3.
	private int maxNumberOfChildren;
	private int numberOfChildren;
	private int childId;  // Number from 0 to 2^dimension
	public String path;  // Sequence of digits showing the path from root to this child
	private Octree parent = null;
	private Octree root = null;
	private Octree[] children = null;
	private OctreeDataItem[] dataItems = null;
	private double[] firstCorner;
	private double[] secondCorner;
	private double[] geometricalCenter;
	private double[] margin;
	private int weight;
	private int leafCutOff = 1;
	private boolean leaf = false;
	private double[] cellCenter;
	private double cellRadius;
	public int levels;
	private double[] size;
	private double[] mid;   // Center of the cell
	private int majorAxis;  // Axis with the maximum spatial extension
	private int minorAxis;  // Axis with the minimal spatial extension
	public static int countChildren = 0;
	public static int countAppBondOp = 0;
	public static int rejectedPaths = 0;


	//
	// Constructors.
	//

	/**
	 *  Construct an octree as a child of "parent" in a tree rooted at
	 *  "root" with given data and corners.
	 */
	public Octree( Octree root, Octree parent, int id, OctreeDataItem[] data,
		double[] firstCorner, double[] secondCorner )
	{
		this.root = root;
		this.childId = id;
		this.path = parent.path + id;
		this.dimension = root.getDimension(); // Assume dimension passed is the same as dimension of data
		// this.dataItems = data;
		this.firstCorner = firstCorner;
		this.secondCorner = secondCorner;
		this.parent = parent;
		this.maxNumberOfChildren = root.getMaxNumberOfChildren( );

		//
		// Set the size in each direction.
		//

		double maxSize = 0.0;
		double minSize = 1.0E6;
		size = new double[ dimension ];
		mid = new double[ dimension ];
		cellRadius = 0.0;
		for ( int j=0; j<dimension; j++ )
		{
			size[j] = Math.abs( secondCorner[j] - firstCorner[j] );
			mid[j] = ( secondCorner[j] + firstCorner[j] ) / 2.0;
			cellRadius += size[j] * size[j];
			if ( size[j] > maxSize )
			{
				maxSize = size[j];
				majorAxis = j;
			}
			if( size[j] < minSize )
			{
				minSize = size[j];
				minorAxis = j;
			}
		}

		cellRadius = Math.sqrt( cellRadius );
		cellRadius /= 2.0;
	}

	/**
	 *  Constructor that is mostly used for instantiating the root tree.
	 */
	public Octree( int spaceDimension, OctreeDataItem[] data, double[] offset )
	{
		//
		// Derive the domain spaned by the coordinates of the elements.
		//

		dimension = spaceDimension; // Assume dimension passed is the same as dimension of data
		dataItems = data;
		weight = data.length;
		margin = new double[ dimension ];

		// System.out.println( "Root tree data length: " + dataItems.length );
		firstCorner = new double[ dimension ];
		secondCorner = new double[ dimension ];
		for ( int i=0; i<dimension; i++ )
		{
			margin[i] = offset[i];
			firstCorner[i] = 1.0E6;
			secondCorner[i] = -1.0E6;
		}

		maxNumberOfChildren = 1;
		for ( int i=0; i<dimension; i++ )
		{
			maxNumberOfChildren <<= 1;
		}

		//
		// Determine the two corners from data coordinates
		//
		
//		System.out.println("data length = " + dataItems.length);

		double x;
		for ( int i=0; i<dataItems.length; i++ )
		{
			if (dataItems[i] == null) continue;
			double coordinate[] = dataItems[i].getCoordinate( );
			for ( int j=0; j<dimension; j++ )
			{
				x = coordinate[j];
				if ( x <= firstCorner[j] ) firstCorner[j] = x;
				if ( x >= secondCorner[j] ) secondCorner[j] = x;
			}
		}

		//
		// Add margins to the corners to make the box somewhat larger than
		// the container of spatial data and set the size in each direction.
		//

		double maxSize = 0.0;
		double minSize = 1.0E6;
		size = new double[dimension];
		cellRadius = 0.0;
		mid = new double[ dimension ];
		for ( int j=0; j<dimension; j++ )
		{
			firstCorner[j] -= margin[j];
			secondCorner[j] += margin[j];
			size[j] = Math.abs( secondCorner[j] - firstCorner[j] );
			mid[j] = ( secondCorner[j] + firstCorner[j] ) / 2.0;
			cellRadius += size[j] * size[j];
			if ( size[j] > maxSize )
			{
				maxSize = size[j];
				majorAxis = j;
			}
			if ( size[j] < minSize )
			{
				minSize = size[j];
				minorAxis = j;
			}
		}

		cellRadius = Math.sqrt( cellRadius );
		cellRadius /=  2.0;

		//
		// We have a root of the octree now
		//

		childId = 0;
		path = "" + childId;
		levels = 0;
		setRoot( this );
		// if ( data.length > 0 ) build( );
	}

	//
	// Methods.
	//

	/**
	 *  Build the whole tree recursively.
	 */
	public void build( )
	{
		countChildren++;


		children = new Octree[ maxNumberOfChildren ];
		Vector[] dataSets = new Vector[ maxNumberOfChildren ]; // Collects data for each child
		// double[] mid = new double[ dimension ];
		double[] childSize = new double[ dimension ];
		double[] corner2;
		int setFirstBit = 1;

		//
		// Determine the center of the parallelepiped.
		//
		
		for ( int i=0; i<dimension; i++ )
		{
			childSize[i] = size[i] / 2.0;
			// mid[i] = ( firstCorner[i] + secondCorner[i] ) / 2;
		}

		//
		// Split the data according to the cell.
		//

		for ( int i=0; i<maxNumberOfChildren; i++ )
		{
			dataSets[i] = new Vector( );
		}

		for ( int i=0; i<dataItems.length; i++ )
		{
			int child = 0;
			if (dataItems[i] == null) continue;
			
			double coordinate[] = dataItems[i].getCoordinate( );
			for ( int j=0; j<dimension; j++ )
			{
				child <<= 1;
				if ( coordinate[j] <= mid[j] )
				{
					child |= setFirstBit;
				}
			}
			dataSets[child].add( dataItems[i] );
		}

		OctreeDataItem[] tmpData = null;
		int dataSize = 0;
		Iterator dataIterator = null;
		int l = 0;
		int testBit;
		double[] secondCorner = null;
		numberOfChildren = 0;

		if ( dataItems.length > leafCutOff )
		{
			for ( int i=0; i<maxNumberOfChildren; i++ )
			{
				dataSize = dataSets[i].size( );
				if ( dataSize > 0 )
				{
					l = 0;
					tmpData = new OctreeDataItem[ dataSize ];
					dataIterator = dataSets[i].iterator( );

					while ( dataIterator.hasNext() )
					{
						tmpData[l] = (OctreeDataItem) dataIterator.next( );
						l++;
					}

					//
					// Calculate the second corner of the child
					//

					corner2 = new double[ dimension ];
					testBit = 1;
					for ( int j=dimension-1; j>=0; j-- )
					{
						if ( (i & testBit) > 0 )
						{
							corner2[j] = mid[j] - childSize[j];
						}
						else
						{
							corner2[j] = mid[j] + childSize[j];
						}
						testBit <<= 1;
					}

					children[i] = new Octree( this.root, this, i, null, mid, corner2 );
					
					children[i].build( tmpData );

				}
				numberOfChildren++;
			}
		}
		else
		{
			setLeafFlag( true );
//			System.out.println( "Warning: Root is also leaf " );
		}
		
		//
		// Set the geometrical center of this node.
		//

		// setGeometricalCenter( );
		// System.out.println( "GC: (" + geometricalCenter[0] + "," + geometricalCenter[1] + "," + geometricalCenter[2] + ")" );
	}

	/**
	 *  Build a subtree recursively.
	 */
	public void build( OctreeDataItem[] data )
	{
		countChildren++;
		weight = data.length;
		


		if ( data.length <= leafCutOff )
		{
			setData( data );
			leaf = true;
			return;
		}

		children = new Octree[ maxNumberOfChildren ];
		Vector[] dataSets = new Vector[ maxNumberOfChildren ]; // Collects data for each child
		// double[] mid = new double[ dimension ];
		double[] childSize = new double[ dimension ];
		double[] corner2;
		int setFirstBit = 1;

		//
		// Determine the center of the parallelepiped.
		//

		for( int i=0; i<dimension; i++ )
		{
			childSize[i] = size[i] / 2.0;
			// mid[i] = ( firstCorner[i] + secondCorner[i] ) / 2;
		}

		//
		// Split the data according to the cell.
		//

		for ( int i=0; i<maxNumberOfChildren; i++ )
		{
			dataSets[i] = new Vector( );
		}

		for ( int i=0; i<data.length; i++ )
		{
			int child = 0;
			double coordinate[] = data[i].getCoordinate( );
			for ( int j=0; j<dimension; j++ )
			{
				child <<= 1;
				if ( coordinate[j] <= mid[j] )
				{
					child |= setFirstBit;
				}
			}
			dataSets[child].add( data[i] );
		}

		OctreeDataItem[] tmpData = null;
		int dataSize = 0;
		Iterator dataIterator = null;
		int l = 0;
		int testBit;
		double[] secondCorner = null;
		numberOfChildren = 0;

		/*
		for( int i=0; i<maxNumberOfChildren; i++ )
		{
			System.out.println( "Data Sizes: " + i + "  " + dataSets[i].size() );
		}
		*/

		for ( int i=0; i<maxNumberOfChildren; i++ )
		{
			dataSize = dataSets[i].size( );
			if ( dataSize > 0 )
			{
				l = 0;
				tmpData = new OctreeDataItem[ dataSize ];
				dataIterator = dataSets[i].iterator( );
				while( dataIterator.hasNext() )
				{
					tmpData[l] = (OctreeDataItem) dataIterator.next( );
					l++;
				}

				//
				// Calculate the second corner of the child.
				//

				corner2 = new double[ dimension ];
				testBit = 1;
				for( int j=dimension-1; j>=0; j-- )
				{
					if ( (i & testBit) > 0 )
					{
						corner2[j] = mid[j] - childSize[j];
					}
					else
					{
						corner2[j] = mid[j] + childSize[j];
					}
					testBit <<= 1;
				}

				children[i] = new Octree( this.root, this, i, null, mid, corner2 );
				children[i].build( tmpData );
				numberOfChildren++;
			}
		}
		

	}

	/**
	 *  Methods to operate with vectors.
	 */
	public final double getDistance( double[] coord1, double[] coord2 )
	{
		//
		// Assume that the two coords have the same dimensionality.
		//

		double distance = 0.0;
		for( int i = 0; i < coord1.length; i++ )
		{
			distance += (coord2[i] - coord1[i]) * (coord2[i] - coord1[i]);
		}
		distance = Math.sqrt( distance );
		return distance;
	}


	/**
	 *  ???
	 */
	public final double[] add( double[] v1, double[] v2 )
	{
		double[] sum = new double[ v1.length ];
		for( int i=0; i<v1.length; i++ )
		{
			sum[i] = v1[i] + v2[i];
		}
		return sum;
	}

	/**
	 *  ???
	 */
	public final double[] scale( double[] v1, double scaleFactor )
	{
		double[] v = new double[ v1.length ];
		for( int i=0; i<v1.length; i++ )
		{
			v[i] = scaleFactor * v1[i];
		}
		return v;
	}

	/**
	 *  ???
	 */
	private void printVector( double[] v )
	{
		System.out.print( "(" );
		for( int i=0; i<v.length-1; i++ )
		{
			System.out.print( " " + v[i] );
			System.out.print( "," );
		}
		System.out.print( " " + v[v.length-1] );
		System.out.println( " )" );
	}

	//
	// Set methods.
	//

	/**
	 *  ???
	 */
	public void setRoot( Octree rootTree )
	{
		this.root = rootTree;
	}

	/**
	 *  ???
	 */
	public void setParent( Octree parentTree )
	{
		this.parent = parentTree;
	}

	/**
	 *  ???
	 */
	public void setData( OctreeDataItem[] data )
	{
		this.dataItems = data;
	}

	/**
	 *  ???
	 */
	public void setLeafFlag( boolean leafFlag )
	{
		this.leaf = leafFlag;
	}

	/**
	 *  ???
	 */
	public void setLeafCutOffCriterion( int leafCutOff )
	{
		this.leafCutOff = leafCutOff;
	}

	/**
	 *  ???
	 */
	public void setChildId( int child )
	{
		childId = child;
	}

	/**
	 *  ???
	 */
	public void setWeight( int weight )
	{
		this.weight = weight;
	}

	//
	// Mixed methods???
	//

	/**
	 *  ???
	 */
	public double[] setGeometricalCenter( )
	{
		geometricalCenter = new double[ dimension ];
		double[] tmp = new double[ dimension ];
		if ( leaf == true )
		{
			for ( int i=0; i<dataItems.length; i++ )
			{
				tmp = add( tmp, dataItems[i].getCoordinate() );
			}
			weight = dataItems.length;
		}
		else
		{
			double[] gc;
			for ( int i=0; i<maxNumberOfChildren; i++ )
			{
				if ( children[i] != null )
				{
					gc = children[i].setGeometricalCenter( );
					tmp = add( tmp, gc );
				}
			}
		}

		for ( int i=0; i<dimension; i++ )
		{
			geometricalCenter[i] = tmp[i];
			geometricalCenter[i] /= weight;
		}

		return tmp;
	}

	//
	// Get methods.
	//

	/**
	 *  ???
	 */
	public Object[] getBonds( double cutOff )
	{
		Vector bondVector = new Vector();
		int firstIndex;

		if ( dataItems != null )
		{
			for ( int i=0; i<dataItems.length; i++ )
			{
				appendBonds( dataItems[i], cutOff, bondVector );
			}
		}
		else
		{
			System.out.println( "This doesn't seem to be a root tree, No data found " );
		}

		bondVector.trimToSize( );
		return bondVector.toArray( );
	}

	//
	// Get methods.
	//

	/**
	 *  Return the vector of all potential covalent bonds.
	 */
	public Vector getBondsVector( double cutOff )
	{
		Vector bondVector = new Vector( );
		double[] firstCoord = new double[ dimension ];
		int firstIndex;

		if ( dataItems != null )
		{
			for( int i=0; i<dataItems.length; i++ )
			{
				appendBonds( dataItems[i], cutOff, bondVector );
			}
		}
		else
		{
			System.out.println( "This doesn't seem to be a root tree, No data found " );
		}
		bondVector.trimToSize( );
		return bondVector;
	}

	/**
	 *  ???
	 */
	public Object[] getHBonds( double cutOff )
	{
		Vector bondVector = new Vector( );
		int firstIndex;

	 	if ( dataItems != null )
		{
			for ( int i=0; i<dataItems.length; i++ )
			{
				appendHBonds( dataItems[i], cutOff, 2, bondVector );
			}
		}
		else
		{
			System.out.println( "This doesn't seem to be a root tree, No data found " );
		}

		bondVector.trimToSize( );
		return bondVector.toArray( );
	}

	/**
	 *  ???
	 */
	public Vector getHBondsVector( double cutOff )
	{
		Vector bondVector = new Vector( );
		int firstIndex;

		if ( dataItems != null )
		{
			for( int i = 0; i < dataItems.length; i++ )
			{
				appendHBonds( dataItems[i], cutOff, 2, bondVector );
			}
		}
		else
		{
			System.out.println( "This doesn't seem to be a root tree, No data found " );
		}

		return bondVector;
	}

    public Vector getHBondInfoVector( double cutOff ) {
	Vector bondList = new Vector();
	double[] firstCoord = new double[dimension];
	int firstIndex;

	//System.out.println ( "Detecting bonds at a cutoff of " + cutOff );
	if( dataItems != null ) {
	    for( int i = 0; i < dataItems.length; i++ ) {
	    	if (dataItems[i] == null) continue;
			firstCoord = dataItems[i].getCoordinate();
			firstIndex = dataItems[i].getIndex();
			appendHBondInfo( firstCoord, firstIndex, cutOff, 2, bondList );
	    }
	}
	else {
	    System.out.println( "This doesn't seem to be a root tree, No data found " );	    
	} 
	return bondList;
    }

    public void appendHBondInfo( double[] coords1, int index1, double cutOff, int indexCutOff, Vector bondList ) {
//	root.countAppBondOp++;
    Octree.countAppBondOp++;
	double distance;
	int index2;
	distance = getDistance( coords1, mid );
	distance -= cellRadius;
	if( distance <= cutOff ) {
	    //System.out.println( "Distance in " + path + "  " + distance );
	    if( leaf ) {
		//System.out.println( "Leaf with " + dataItems.length + " components" );
		for( int j = 0; j < dataItems.length; j++ ) {
			if (dataItems[j] == null) continue;
			
		    if( getDistance( dataItems[j].getCoordinate(), coords1 ) <= cutOff ) {
			//System.out.println( "Appending bonds in " + path );
			index2 = dataItems[j].getIndex();
			if( index1 < index2 - indexCutOff ) {
			    bondList.add( new BondInfo( index1, dataItems[j].getIndex() ) );
			}
		    }
		}
	    }
	    else {
		for( int child = 0; child < maxNumberOfChildren; child++ ) {
		    if( children[child] != null ) {
			children[child].appendHBondInfo( coords1, index1, cutOff, indexCutOff, bondList );
		    }
		}
	    }
	}
	else {
	    //root.rejectedPaths++;
	    //System.out.println( "Rejecting Path: " + path );
	    return;
	}
	return;
    }
	
        /**
	 *  ???
	 */
	public void appendBonds( OctreeDataItem dataItem, double cutOff,
		Vector bondVector )
	{
		Octree.countAppBondOp++;
		
		double distance;
		int index1, index2;

		double[] coords1 = dataItem.getCoordinate();
		index1 = dataItem.getIndex();
		distance = getDistance( coords1, mid );
		distance -= cellRadius;

		if ( distance > cutOff ) return;

		if ( leaf )
		{
			for ( int j=0; j<dataItems.length; j++ )
			{
				if ( getDistance( dataItems[j].getCoordinate(), coords1 ) <= cutOff )
				{
					index2 = dataItems[j].getIndex( );
					if ( index1 < index2 )
					{
						bondVector.add( new Bond( ((OctreeAtomItem)dataItem).getAtom(), ((OctreeAtomItem)dataItems[j]).getAtom() ) );
					}
				}
			}
		}
		else
		{
			for ( int child=0; child<maxNumberOfChildren; child++ )
			{
				if ( children[child] != null )
				{
					children[child].appendBonds( dataItem, cutOff, bondVector );
				}
			}
		}
	}

	/**
	 *  ???
	 */
	public void appendHBonds( OctreeDataItem dataItem, double cutOff,
		int indexCutOff, Vector bondVector )
	{
		Octree.countAppBondOp++;
		
		double distance;
		int index1, index2;

		double[] coords1 = dataItem.getCoordinate();
		index1 = dataItem.getIndex();
		distance = getDistance( coords1, mid );
		distance -= cellRadius;

		if ( distance > cutOff ) return;

		if ( leaf )
		{
			for ( int j=0; j<dataItems.length; j++ )
			{
				if ( getDistance( dataItems[j].getCoordinate(), coords1 ) <= cutOff )
				{
					index2 = dataItems[j].getIndex();
					if ( index1 < index2 - indexCutOff )
					{
						bondVector.add( new Bond( ((OctreeAtomItem)dataItem).getAtom(), ((OctreeAtomItem)dataItems[j]).getAtom() ) );
					}
				}
			}
		}
		else
		{
			for ( int child=0; child<maxNumberOfChildren; child++ )
			{
				if ( children[child] != null )
				{
					children[child].appendHBonds( dataItem, cutOff, indexCutOff, bondVector );
				}
			}
		}
	}

	/**
	 *  ???
	 */
	public  int getDimension( )
	{
		return dimension;
	}

	/**
	 *  ???
	 */
	public int getMaxNumberOfChildren( )
	{
		return maxNumberOfChildren;
	}

	/**
	 *  ???
	 */
	public int getNumberOfChildren( )
	{
		return numberOfChildren;
	}

	/**
	 *  ???
	 */
	public int getWeight( )
	{
		return weight;
	}

	/**
	 *  ???
	 */
	public void printChildNamingScheme( )
	{
	}
}