#ifndef HASH_INTERFACE_H
#define HASH_INTERFACE_H

#include "GeomHash.h"
#include "RigidTrans3.h"
#include "Interface.h"

#define defaultThreshold 4.5
#define hashParam_d 3
#define hashParam_cs 0.5

/*
CLASS
  HashInterface

  The class can build Interface object efficiently by using Geometric Hashing
  for the particles of the first molecule.
  
KEYWORDS
  HashInterface, Interface, GeomHash, RigidTrans3

AUTHORS
  Inbal Halperin (inbalhal@math.tau.ac.il), 
  Dina Duhovny (duhovka@math.tau.ac.il) and 
  Yuval Inbar (inbaryuv@math.tau.ac.il)

  Copyright: SAMBA group, Tel-Aviv Univ. Israel, 1997-9.

CHANGES LOG

<UL>
  <LI>
23.11.03 Dina
Implementation of probe interface definition.
  <\LI>
  <LI>
23.11.03 Dina
Option to compute together two Interface objects (one for each molecule)
 <\LI>
</UL>

GOALS
  Create an Interface for two Molecule descendants using Geometric
  Hashing. This implementation allows fast interface generation.

  The class supports 2 interface definitions:
  1. If the distance between two atoms, one from each of the molecules, is less
  than distance_threshold, the two corresponding atoms are marked as interface
  atom pair.
  2. If the distance between two atoms, one from each of the molecules, is less
  than the sum of their vdW radii and probe diameter, the two corresponding
  atoms are marked as interface atom pair. Intutively, two atoms are interface
  pair, if the probe sphere can't be placed between them.
  
USAGE
  In the construction the first molecule is preprocessed into GeomHash. This
  allows fast mapping of interface pairs.

  The class has several options to build the interface. First the user should
  pick the interface definition he needs (see Goals). For each of the two
  definitions, 8 functions with different arguments are provided. Below we
  summarize different options:

  The first molecule remains the same in all cases. The second molecule can be
  the same that was given in constructor or can be specified by the user. In
  some cases, like generating docking results interfaces, the second molecule
  is always the same, however the orientation is different. In this case the
  user may specufy the 3D transformation to apply to the second molecule.

  Since the Interface object is not symmetric, it is possible to specify which
  molecule will be used as first index in the Interface. The first molecule
  MOL1 is the default. Also it is possible to generate both interfaces at once.

  The following example demonstrates how two Interfaces are  created for 
  one complex. The complex in the example is composed of an Antibody=Ab and an
  Antigen=Ag. Ab and Ab are of type Molecule<Atom>. 
  
  EXAMPLE  
  HashInterface<Molecule<Atom> > hashInterface(Ab,Ag);
  Interface face1, face2;
  hashInterface.buildInterface(face1, face2);
  for(Interface::iterator it=face2.begin(); it!=face2.end(); ++it) {
    cout << "serial:  " << *it << endl;    
    cout << "pdb first atom:  " << Ab[(*it).second()]<<endl;
    cout << "pdb second atom:  " << Ag[(*it).first()]<<endl;
  }
  END
*/

template< class MoleculeT>
class HashInterface {
 public:
  
  enum MOL_INDEX{ MOL1, MOL2 };

  // GROUP: Constructors

  //// Construct hash interface, hash first molecule
  HashInterface(const MoleculeT& iMol1, const MoleculeT& iMol2,
                const unsigned short d=hashParam_d, const float cs=hashParam_cs);

  // GROUP: distance defined interface
  
  //// Construct interface for two molecules.
  void buildInterface(Interface& face,
		      const float distThr=defaultThreshold, const MOL_INDEX molIndex=MOL1);

  //// Construct interface for two molecules. Apply transformation on a second
  //   molecule.
  void buildInterface(Interface& face, const RigidTrans3& trans2,
		      const float distThr=defaultThreshold, const MOL_INDEX molIndex=MOL1);

  //// Construct interface for two molecules: first from class object and the
  //   second given as a parameter.
  void buildInterface(Interface& face, const MoleculeT& iMol2,
		      const float distThr=defaultThreshold, const MOL_INDEX molIndex=MOL1);
 
  //// Construct interface for two molecules: first from class object and the
  //   second given as a parameter. Apply transformation on a second
  //   molecule.
  void buildInterface(Interface& face, const MoleculeT& iMol2, const RigidTrans3& trans2,
		      const float distThr=defaultThreshold, const MOL_INDEX molIndex=MOL1);

  //// Construct two interfaces for two molecules.
  void buildInterface(Interface& face1, Interface& face2, const float distThr=defaultThreshold);

  //// Construct two interfaces for two molecules. Apply transformation on a second
  //   molecule.
  void buildInterface(Interface& face1, Interface& face2,
		      const RigidTrans3& trans2, const float distThr=defaultThreshold);

  //// Construct two interfaces for two molecules: first from class object and the
  //   second given as a parameter.
  void buildInterface(Interface& face1, Interface& face2,
		      const MoleculeT& iMol2, const float distThr=defaultThreshold);
 
  //// Construct two interfaces for two molecules: first from class object and
  //   the second given as a parameter. Apply transformation on a second
  //   molecule.  
  void buildInterface(Interface& face1, Interface& face2, const MoleculeT& iMol2,
		      const RigidTrans3& trans2, const float distThr=defaultThreshold);

  // GROUP: probe sphere defined interface
  
  //// Construct interface for two molecules.
  void buildProbeInterface(Interface& face, const float probeRadius, const MOL_INDEX molIndex=MOL1);

  //// Construct interface for two molecules. Apply transformation on a second
  //   molecule.
  void buildProbeInterface(Interface& face, const RigidTrans3& trans2,
			   const float probeRadius, const MOL_INDEX molIndex=MOL1);

  //// Construct interface for two molecules: first from class object and the
  //   second given as a parameter.
  void buildProbeInterface(Interface& face, const MoleculeT& iMol2,
			   const float probeRadius, const MOL_INDEX molIndex=MOL1);
 
  //// Construct interface for two molecules: first from class object and the
  //   second given as a parameter. Apply transformation on a second
  //   molecule.
  void buildProbeInterface(Interface& face, const MoleculeT& iMol2, const RigidTrans3& trans2,
			   const float probeRadius, const MOL_INDEX molIndex=MOL1);

  //// Construct two interfaces for two molecules.
  void buildProbeInterface(Interface& face1, Interface& face2, const float probeRadius);

  //// Construct two interfaces for two molecules. Apply transformation on a second
  //   molecule.
  void buildProbeInterface(Interface& face1, Interface& face2,
			   const RigidTrans3& trans2, const float probeRadius);

  //// Construct two interfaces for two molecules: first from class object and the
  //   second given as a parameter.
  void buildProbeInterface(Interface& face1, Interface& face2,
			   const MoleculeT& iMol2, const float probeRadius);
 
  //// Construct two interfaces for two molecules: first from class object and
  //   the second given as a parameter. Apply transformation on a second
  //   molecule.  
  void buildProbeInterface(Interface& face1, Interface& face2, const MoleculeT& iMol2,
			   const RigidTrans3& trans2, const float probeRadius);

 protected:
  //// function that computes interface pairs and stores them in the Interface
  //   object using distance threshold definition
  void internalBuildInterface(Interface& face, const MoleculeT& secondMol,
			      const float distThr=defaultThreshold, const MOL_INDEX molIndex=MOL1);

  //// function that computes interface pairs and stores them in two Interface
  //   objects (molecule1 and molecule2) using distance threshold definition
  void internalBuildInterface(Interface& face1, Interface& face2,
			      const MoleculeT& secondMol, const float distThr=defaultThreshold);
  
  //// function that computes interface pairs and stores them in the Interface
  //   object using probe sphere definition
  void internalBuildProbeInterface(Interface& face, const MoleculeT& secondMol,
				   const float probeRadius, const MOL_INDEX molIndex=MOL1);

  //// function that computes interface pairs and stores them in the Interface
  //   object using probe sphere definition
  void internalBuildProbeInterface(Interface& face, Interface& face2,
				   const MoleculeT& secondMol, const float probeRadius);
  
 protected:
  //// First molecule
  const MoleculeT& mol1;

  //// Second molecule
  const MoleculeT& mol2;

  //// GeomHash for fast atomic search
  GeomHash< Vector3, unsigned int> gHash1;
};

template< class MoleculeT>
HashInterface< MoleculeT>::HashInterface(const MoleculeT& iMol1,	const MoleculeT& iMol2,
					const unsigned short d, const float cs) 
  : mol1(iMol1), mol2(iMol2), gHash1(d,cs) {
  for(unsigned int i=0; i<mol1.size(); i++) 
    gHash1.insert(mol1[i].position(), i);
}

/*************************** Probe interface functions ********************/

template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face, const float probeRadius,
						   const MOL_INDEX molIndex) {
  internalBuildProbeInterface(face, mol2, probeRadius, molIndex);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face, const RigidTrans3& trans2,
						   const float probeRadius, const MOL_INDEX molIndex) {
  MoleculeT transedMol2=mol2;
  transedMol2.rigidTrans(trans2);
  internalBuildProbeInterface(face, transedMol2, probeRadius, molIndex);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face, const MoleculeT& iMol2,
						   const float probeRadius, const MOL_INDEX molIndex) {
  internalBuildProbeInterface(face, iMol2, probeRadius, molIndex);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face, const MoleculeT& iMol2,
						   const RigidTrans3& trans2,
						   const float probeRadius, const MOL_INDEX molIndex) {
  MoleculeT transedMol2=iMol2;
  transedMol2.rigidTrans(trans2);
  internalBuildProbeInterface(face, transedMol2, probeRadius, molIndex);
}
  
template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face1, Interface& face2, const float probeRadius) {
  internalBuildProbeInterface(face1, face2, mol2, probeRadius);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face1, Interface& face2,
						   const MoleculeT& iMol2, const float probeRadius) {
  internalBuildProbeInterface(face1, face2, iMol2, probeRadius);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face1, Interface& face2,
						   const RigidTrans3& trans2, const float probeRadius) {
  MoleculeT transedMol2=mol2;
  transedMol2.rigidTrans(trans2);
  internalBuildProbeInterface(face1, face2, transedMol2, probeRadius);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildProbeInterface(Interface& face1, Interface& face2, const MoleculeT& iMol2,
						   const RigidTrans3& trans2, const float probeRadius) {
  MoleculeT transedMol2=iMol2;
  transedMol2.rigidTrans(trans2);
  internalBuildProbeInterface(face1, face2, transedMol2, probeRadius);
}

/*************************** Distance interface functions ********************/

template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face, const float distThr, const MOL_INDEX molIndex) {
  internalBuildInterface(face, mol2, distThr, molIndex);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face, const RigidTrans3& trans2,
					      const float distThr, const MOL_INDEX molIndex) {
  MoleculeT transedMol2=mol2;
  transedMol2.rigidTrans(trans2);
  internalBuildInterface(face, transedMol2, distThr, molIndex);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face, const MoleculeT& iMol2,
					      const float distThr, const MOL_INDEX molIndex) {
  internalBuildInterface(face, iMol2, distThr, molIndex);
}


template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face, const MoleculeT& iMol2, const RigidTrans3& trans2,
					      const float distThr, const MOL_INDEX molIndex) {
  MoleculeT transedMol2=iMol2;
  transedMol2.rigidTrans(trans2);
  internalBuildInterface(face, transedMol2, distThr, molIndex);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face1, Interface& face2, const float distThr) {
  internalBuildInterface(face1, face2, mol2, distThr);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face1, Interface& face2,
					      const RigidTrans3& trans2, const float distThr) {
  MoleculeT transedMol2=mol2;
  transedMol2.rigidTrans(trans2);
  internalBuildInterface(face1, face2, transedMol2, distThr);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face1, Interface& face2,
					      const MoleculeT& iMol2, const float distThr) {
  internalBuildInterface(face1, face2, iMol2, distThr);
}

template< class MoleculeT>
void HashInterface< MoleculeT>::buildInterface(Interface& face1, Interface& face2, const MoleculeT& iMol2,
					      const RigidTrans3& trans2, const float distThr) {
  MoleculeT transedMol2=iMol2;
  transedMol2.rigidTrans(trans2);
  internalBuildInterface(face1, face2, transedMol2, distThr);
}

/*************************** Build interface functions ********************/

template< class MoleculeT>
void HashInterface< MoleculeT>::internalBuildInterface(Interface& face, const MoleculeT& secondMol,
						      const float distThr, const MOL_INDEX molIndex) {
  float distThr2 = distThr*distThr;
  for(unsigned int secondMolIndex=0; secondMolIndex < secondMol.size(); secondMolIndex++) {
    HashResult<unsigned int> hashResult;
    gHash1.query(secondMol(secondMolIndex), distThr, hashResult);
    for(HashResult<unsigned int>::iterator targetIter = hashResult.begin();
	targetIter!=hashResult.end(); targetIter++) {
      float distance2 = mol1(*targetIter).dist2(secondMol(secondMolIndex));
      if(distance2 < distThr2) {
	if(molIndex==MOL1)
	  face.addAdjacency(*targetIter, secondMolIndex, sqrt(distance2));
	else
	  face.addAdjacency(secondMolIndex, *targetIter, sqrt(distance2));
      }
    }
  }
}

template< class MoleculeT>
void HashInterface< MoleculeT>::internalBuildInterface(Interface& face1, Interface& face2,
						      const MoleculeT& secondMol, const float distThr) {
  float distThr2 = distThr*distThr;
  for(unsigned int secondMolIndex=0; secondMolIndex < secondMol.size(); secondMolIndex++) {
    HashResult<unsigned int> hashResult;
    gHash1.query(secondMol(secondMolIndex), distThr, hashResult);
    for(HashResult<unsigned int>::iterator targetIter = hashResult.begin();
	targetIter!=hashResult.end(); targetIter++) {
      float distance2 = mol1(*targetIter).dist2(secondMol(secondMolIndex));
      if(distance2 < distThr2) {
	face1.addAdjacency(*targetIter, secondMolIndex, sqrt(distance2));
	face2.addAdjacency(secondMolIndex, *targetIter, sqrt(distance2));
      }
    }
  }
}

template< class MoleculeT>
void HashInterface< MoleculeT>::internalBuildProbeInterface(Interface& face1, Interface& face2,
							   const MoleculeT& secondMol, const float probeRadius) {
  float probeDiameter = 2*probeRadius;
  for(unsigned int secondMolIndex=0; secondMolIndex < secondMol.size(); secondMolIndex++) {
    float radius = secondMol[secondMolIndex].getRadius() + probeDiameter + secondMol[secondMolIndex].maxRadius;
    HashResult<unsigned int> hashResult;
    gHash1.query(secondMol(secondMolIndex), radius, hashResult);
    for(HashResult<unsigned int>::iterator targetIter = hashResult.begin();
	targetIter!=hashResult.end(); targetIter++) {
      float distance2 = mol1(*targetIter).dist2(secondMol(secondMolIndex));
      float maxDist = secondMol[secondMolIndex].getRadius() + mol1[*targetIter].getRadius() + probeDiameter;
      if(distance2 < maxDist*maxDist) {
	face1.addAdjacency(*targetIter, secondMolIndex, sqrt(distance2));
	face2.addAdjacency(secondMolIndex, *targetIter, sqrt(distance2));
      }
    }
  }
}
  
template< class MoleculeT>
void HashInterface< MoleculeT>::internalBuildProbeInterface(Interface& face, const MoleculeT& secondMol,
							   const float probeRadius, const MOL_INDEX molIndex) {
  float probeDiameter = 2*probeRadius;
  for(unsigned int secondMolIndex=0; secondMolIndex < secondMol.size(); secondMolIndex++) {
    float radius = secondMol[secondMolIndex].getRadius() + probeDiameter + secondMol[secondMolIndex].maxRadius;
    HashResult<unsigned int> hashResult;
    gHash1.query(secondMol(secondMolIndex), radius, hashResult);
    for(HashResult<unsigned int>::iterator targetIter = hashResult.begin();
	targetIter!=hashResult.end(); targetIter++) {
      float distance2 = mol1(*targetIter).dist2(secondMol(secondMolIndex));
      float maxDist = secondMol[secondMolIndex].getRadius() + mol1[*targetIter].getRadius() + probeDiameter;
      if(distance2 < maxDist*maxDist) {
	if(molIndex==MOL1)
	  face.addAdjacency(*targetIter, secondMolIndex, sqrt(distance2));
	else
	  face.addAdjacency(secondMolIndex, *targetIter, sqrt(distance2));
      }
    }
  }
}

#endif