// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//

#ifndef BALL_STRUCTURE_BINARYFINGERPRINTMETHODS_H
#define BALL_STRUCTURE_BINARYFINGERPRINTMETHODS_H

#ifndef BALL_COMMON_H
#       include <BALL/common.h>
#endif

#ifndef BALL_DATATYPE_OPTIONS_H
#       include <BALL/DATATYPE/options.h>
#endif


#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/incremental_components.hpp>
#include <boost/thread/mutex.hpp>
#include <boost/thread/thread.hpp>
#include <boost/unordered_set.hpp>


#include <map>
#include <set>
#include <vector>


namespace BALL 
{
	/** 	Binary Fingerprint Methods
		\brief This class provides efficient similarity calculation functionality for 2D binary fingerprints.
		
		BinaryFingerprintMethods uses a blocked inverted index data structures as a representation of 2D binary fingerprints.
		This data structure enables efficient calculation of similarity coefficients which use shared feature counts
		between two fingerprints. The shared feature count is the number of on-bits which two different fingerprints have
		in common. Currently the class provides the Tanimoto coefficient. \n \n
		Based on this similarity algorithm the class provides different useful public methods: \n
		\link BinaryFingerprintMethods::cutoffSearch Cutoff Search \endlink: \n
		Similarities of all molecules in a query library to all molecules in a target library are calculated and the pairs whose similarities
		exceed a predefined cutoff are written to an output file. The performance of the method depends on size of the input query. For a single
		molecule query, the method has minimum performance and grows exponentially with the number of query molecules reaching a maximum performance
		at about 400 query molecules. As a consequence it is faster to bundle cutoff searches for multiple molecules if available. \n
		\link BinaryFingerprintMethods::connectedComponents Connected Components \endlink: \n
		For an input library all pairwise similarities are calculated. A similarity graph is generated by inserting molecule pairs as edges.
		A variable similarity cutoff is applied to insert only those edges whose corresponding similarity exceeds this cutoff. Thereby, a
		similarity network is created and the connected components are finally calculated. If desired, the nearest neighbour of every molecule
		can additionally be returned. \n
		\link BinaryFingerprintMethods::averageLinkageClustering Hierarchical Clustering \endlink: \n
		A set of input molecules is hierarhically clustered according to the average linkage criterion. Two different algorithms are implemented
		to make efficient use of the inverted index algorithm and the available hardware resources. Both methods are based on %Reciprocal Nearest
		Neighbours. The \link BinaryFingerprintMethods::averageLinkageParallel parallel variant \endlink is used to handle large input library sizes. 
		If the similarity matrix for the remaining clusters or possibly all input molecules fit into main memory, the matrix is calculated and 
		the \link BinaryFingerprintMethods::NNChainCore Nearest Neighbour Chain algorithm \endlink is used. \n 
		Finally, the method of \link BinaryFingerprintMethods::clusterSelectionKGS Kelly et al. \endlink is used to select an level for cluster selection 
		and the created clusters are returned. \n
		\link BinaryFingerprintMethods::calculateSelectionMedoid Calculate Medoid \endlink : \n
		For a set of input molecules the Medoid is calculated. The medoid is the molecule with the highest averaged similarity to all other molecules.
		Additionally, the method returns the averaged pairwise similarities for every molecule in the input set. \n \n
		Citations:\n
		Clustering Algorithms: F. Murtagh, Compstat Lectures vol. 4, 1985, Physica-Verlag Würzburg-Wien. (Volume titel: Multidimensional clustering algorithms). \n
		Similarity Update: G. Lance and W. Williams, Comput J. (1967), 9, 373-380. (doi: 10.1093/comjnl/9.4.373).\n
		Cutoff Selection: L.A. Kelley, S.P. Gardner and M.J. Sutcliffe, %Protein Eng. (1996), 9, 1063-1065. (doi: 10.1093/protein/9.11.1063).
	 */
	class BALL_EXPORT BinaryFingerprintMethods
	{
		struct InvertedIndex;
		
		/**
		 * @name Constant Definitions
		 */
		//@{
		
		typedef std::vector<InvertedIndex*> InvertedIndices;
		typedef std::vector<std::vector<unsigned short> > FingerprintFeatures;
		
		typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> SimilarityGraph;
		typedef boost::graph_traits<SimilarityGraph>::vertex_descriptor Vertex;
		typedef boost::graph_traits<SimilarityGraph>::vertices_size_type VertexIndex;
		typedef boost::disjoint_sets<Vertex*, VertexIndex*, boost::find_with_full_path_compression> DisjointSet;
		
		//@}
		
		public:
			/**
			 * @name Constant Definitions
			 */
			//@{
			
			/**
			 * Options
			 */
			struct BALL_EXPORT Option
			{
				/**
				 * The number of molecules which are grouped together to form a single inverted index data structure.
				 */
				static const String BLOCKSIZE;
				
				/**
				 * The similarity cutoff which is used for connected components decomposition or similarity search.
				 */
				static const String SIM_CUTOFF;
				
				/**
				 * Floating point precision for comparisons of similarity coefficents.
				 */
				static const String PRECISION;
				
				/**
				 * Boolean flag indicates if nearest neighbour information is stored during connected components decomposition.
				 */
				static const String STORE_NN;
				
				/**
				 * Maximal number of clusters for which the similarity half-matrix is calculated in the clustering step.
				 * This value is by default set to 10000 and larger values have to be HANDLED WITH CARE due to quadratic memory requirements.
				 */
				static const String MAX_CLUSTERS;
				
				/**
				 * Number of parallel threads to be used. Parallelization via BOOST threads.
				 */
				static const String N_THREADS;
				
				/**
				 * Intensity level of verbose output
				 */
				static const String VERBOSITY;
			};
			
			/**
			 * Default values
			 */
			struct BALL_EXPORT Default
			{
				static const unsigned short BLOCKSIZE;
				static const float SIM_CUTOFF;
				static const float PRECISION;
				static const unsigned int MAX_CLUSTERS;
				static const bool STORE_NN;
				static const unsigned int N_THREADS;
				static const int VERBOSITY;
			};
			
			
			/**
			* @name Constructors and Destructors 
			*/
			//@{
			
			/**
			* Default constructor.
			*/
			BinaryFingerprintMethods();
			
			
			/**
			* Constructor
			* @param options Options to use. Defaults will be overwritten. 
			*/
			BinaryFingerprintMethods(const Options& options);
			
			
			/**
			 * Constructor
			 * @param options Options to use. Defaults will be overwritten.
			 * @param lib_features Target library as a lists of integer fingerprint features.
			 */
			BinaryFingerprintMethods(const Options& options, const FingerprintFeatures& lib_features);
			
			
			/**
			 * Constructor
			 * @param options Options to use. Defaults will be overwritten.
			 * @param lib_features Target library as lists of integer fingerprint features.
			 * @param query_features Query library as lists of integer fingerprint features.
			 */
			BinaryFingerprintMethods(const Options& options, const FingerprintFeatures& lib_features, const FingerprintFeatures& query_features);
			
			
			/**
			 * Copy constructor.
			 */
			BinaryFingerprintMethods(const BinaryFingerprintMethods& bfm);
			
			
			/**
			* Destructor.
			*/
			virtual ~BinaryFingerprintMethods();
			
			//@}
			
			/** 
			 * @name Assignment 
			 */
			//@{
			
			/** 
			 * Assignment operator
			 */
			BinaryFingerprintMethods& operator = (const BinaryFingerprintMethods& bfm);
			
			//@}
			
			
			/** 
			 * Wrapper for different binary fingerprint parsers which returns a vector of integer features.
			 * The returned features lie in the interval [1, max_feature_id + 1]. Thus, no feature has ID = 0 which is important for inverted index implementation.
			 * @param fprint The fingerprint as a separated list of integer features.
			 * @param features A vector reference to return the fingerprint features. The vector will be cleaned first.
			 * @param fp_type Fingerprint encoding: 1 = Binary Bitstring, 2 = Separated list of integer features
			 * @param delim The delimiter character of the list.
			 * @return True if feature list has been parsed successful.
			 */
			static bool parseBinaryFingerprint(const String& fprint, std::vector<unsigned short>& features, unsigned int fp_type, const char* delim=",");
			
			
			/** 
			 * @name Accessors
			 */
			//@{
			
			/**
			 * Add a target library containing binary fingerprints as feature lists.
			 * @param lib_features Target library as lists of integer fingerprint features.
			 * @return True if setup was successful, false if formats of input features are invalid.
			 */
			bool setLibraryFeatures(const FingerprintFeatures& lib_features);
			
			
			/**
			 * Add a query library containing binary fingerprints as feature lists.
			 * @param query_features Query library as lists of integer fingerprint features.
			 * @return True if setup was successful, false if formats of input features are invalid.
			 */
			bool setQueryFeatures(const FingerprintFeatures& query_features);
			
			
			/**
			 * Get number of molecules in target library.
			 * @return Number of target molecules.
			 */
			unsigned int getTargetLibrarySize() const;
			
			
			/**
			 * Get number of molecules in query library.
			 * @return Number of query molecules.
			 */
			unsigned int getQueryLibrarySize() const;
			
			
			/**
			 * Get const reference to options.
			 * @return Const reference to options of BinaryFingerprintMethods.
			 */
			const Options& getOptions() const;
			
			
			/** 
			 * Set verbosity level.
			 * @param verbosity verbosity level to use.
			 */
			void setVerbosityLevel(const int verbosity);
			
			//@}
			
			
			/** 
			 * @name Applications
			 */
			//@{
			
			/**
			 * Perform similarity search for molecules of a query library in a target library and keep pairs with similarity exceeding sim_cutoff.
			 * @param sim_cutoff Similarity cutoff to apply for cutoff search, i.e. only keep similarities above cutoff.
			 * @param file_name Output file name to store result pairs with similarities as space separated CSV file.
			 */
			bool cutoffSearch(const float sim_cutoff, const String& outfile_name);
			
			
			/** 
			 * Calculate similarity network of a set of molecules defined by selection and return the connected components.
			 * @param selection Indices of molecules in lib_features_ array for which connected components analysis should be performed.
			 * @param ccs Data structure to return connected components. Eeach ccs member is a single connected component which stores the indices
			 * (positions in the vector) of the connected component members in selection. 
			 * @param nn_data If store_nns = true, nn_data will store the nearest neighbour information (index into selection and similarity)
			 * for all connected component members.
			 * @param sim_cutoff Similarity cutoff to apply for similarity network generation, i.e. only keep similarities above cutoff.
			 * @param store_nns If true, the nearest neighbour and its similarity are stored during connected components calculation.
			 */
			bool connectedComponents(const std::vector<unsigned int>& selection, 
						 std::vector<std::vector<unsigned int> >& ccs, 
						 std::vector<std::vector<std::pair<unsigned int, float> > >& nn_data, 
						 const float cutoff,
						 const bool store_nns = false);
			
			
			/** 
			 * Average linkage clustering of a set of molecules defined by selection.
			 * @param selection Indices of molecules in lib_features_ vector for which connected components analysis should be performed
			 * @param nn_data Nearest neighbour information for molecules in selection if already calculated, otherwise an empty vector.
			 * @param cluster_selection Final clustering, i.e. a cluster ID is assigned to every leaf node.
			 * @return True if clustering was successful, false otherwise.
			 */
			bool averageLinkageClustering(const std::vector<unsigned int>& selection, 
						      std::vector<std::pair<unsigned int, float> >& nn_data, 
						      std::map<unsigned int, std::vector<unsigned int> >& cluster_selection);
			
			
			/** 
			 * Calculation of the medoid of a set of molecules defined by selection. Medoid has the highest average similarity to all other compounds.
			 * @param selection Indices of molecules in lib_features_ vector for which the medoid should be calculated.
			 * @param medoid_index Reference to return the index of the calculated medoid in the passed molecule indices (selection).
			 * @param avg_sims Reference to vector which finally will store the average similarity of every molecule to all others in selection.
			 * @return True if calculation has been successful, otherwise false.
			 */
			bool calculateSelectionMedoid(const std::vector<unsigned int>& selection, unsigned int& medoid_index, std::vector<float>& avg_sims);
			
			//@}
			
		private:
			/**
			 * Struct which stores information for a single BOOST thread.
			 */
			struct ThreadData
			{
				LongSize first;
				LongSize last;
				
				float nn_sim;
				unsigned int blocksize;
				unsigned int dataset_size;
				unsigned int active_iids_size;
				
				unsigned short* cc_row;
				unsigned short** cc_matrix;
				
				float* float_array;
				double** dprec_matrix;
				unsigned int* uint_array;
				
				String outfile_name;
			};
			
			
			/**
			 * Core data structure of an InvertedIndex.
			 * For every unique feature in an inverted index, a FeatureList is created.
			 * The FeatureList contains an integer array which stores the position of every molecule in the corresponding InvertedIndex which shares this feature.
			 */ 
			struct FeatureList
			{
				// FeatureID which corresponds in principle to the index of a bit in a 2D fingerprint.
				unsigned short feature_id;
				
				// Array which stores the sorted list of all molecules which all share the associated feature_id.
				// Molecules are only represented by their positions in the Block they belong to.
				unsigned short* block_positions;
			};
			
			
			/**
			 * InvertedIndex core data structure.
			 */
			struct InvertedIndex
			{
				// Number of molecules in the block.
				unsigned int n_molecules;
				
				// Array which stores the number of features of every molecule in the InvertedIndex
				unsigned short* n_features;
				
				// Array which stores the ID of the parent cluster to which the molecule at the InvertedIndex position belongs to.
				unsigned int* parent_clusters;
				
				// Array of FeatureLists implemented as a skip list.
				FeatureList* feature_skip_list;
			};
			
			
			/**
			 * Struct which stores information of a single cluster.
			 */
			struct Cluster
			{
				// Molecule ID for leaf nodes. Otherwise internal Cluster ID.
				unsigned int c_id;
				
				// Cluster index [0, n_molecules-1] is used to map from 2D matrix to 1D array.
				unsigned int c_index;
				
				// Number of leaf nodes inside in the corresponding cluster.
				unsigned int c_size;
				
				// Left child. NULL if cluster is a leaf.
				Cluster* l_child;
				
				// Right child. NULL if cluster is a leaf.
				Cluster* r_child;
				
				// Pointer to the clusters current nearest neigbour.
				Cluster* nn;
				
				// Flag for various purposes.
				bool is_rnn;
				
				// Sum of all intra-cluster pairwise similarities within this cluster. Used finally for KGS penalty calculation.
				double sim_sum;
				
				// Pointer to predecessor cluster in the Nearest Neighbour Chain. NULL if not member of the Nearest Neighbour Chain.
				Cluster* predecessor;
				
				// Similarity to predecessor in Nearest Neighbour Chain. 0.0 if not member of Nearest Neighbour Chain.
				float predecessor_sim;
				
				// Sum of all inter-cluster pairwise similarities of clusters nn_chain[i] and nn_chain[i-1]. 0.0 if not member of the NNChain.
				double predecessor_sim_sum;
				
				// Inverted indices of all leaf molecules within this cluster.
				boost::unordered_set<InvertedIndex*> c_members;
				
				// Vector of pointers to FingerprintFeatures in lib_features_ of all leaf molecules within this cluster.
				std::vector<const std::vector<unsigned short>* > leaf_members;
			};
			
			
			enum clustering_methods_
			{
				STORED_DATA_PARALLEL,
				STORED_MATRIX
			};
			
			
			/**
			 * BOOST DisjointSet used for incremental connected components calculation.
			 */
			DisjointSet* ds;
			
			
			/**
			 * BOOST parent data structure used for incremental connected components calculation.
			 */
			std::vector<Vertex> parent;
			
			
			/**
			 * BOOST rank data structure used for incremental connected components calculation.
			 */
			std::vector<VertexIndex> rank;
			
			
			/**
			 * General options for BinarySimilarityMethods.
			 */
			Options options_;
			
			
			/**
			 * Number of parallel BOOST threads to use.
			 */
			unsigned int n_threads_;
			
			
			/**
			 * Target input library stored as a vector of feature vectors.
			 * A single feature vector is a strictly decreasing list of fingerprint feature IDs (> 0) and a terminating 0.
			 */
			const FingerprintFeatures* lib_features_;
			
			
			/**
			 * Vector of pointers to InvertedIndices for library molecules.
			 */
			InvertedIndices lib_iindices_;
			
			
			/**
			 * Query input library stored as a vector of feature vectors
			 * A single feature vector is a strictly decreasing list of fingerprint feature IDs (> 0) and a terminating 0.
			 */
			const FingerprintFeatures* query_features_;
			
			
			/**
			 * Vector of pointers to InvertedIndices for query molecules.
			 */
			InvertedIndices query_iindices_;
			
			
			/**
			 * Array to store executed BOOST threads.
			 */
			boost::thread* threads_;
			
			
			/**
			 * Boost mutual exclusion variable to secure parallel execution during different tasks.
			 */
			boost::mutex out_mutex_;
			
			
			/**
			 * Array to store ThreadData for every thread.
			 */
			ThreadData* thread_data_;
			
			
			/**
			 * Size of InvertedIndices. I.e. the numbers of molecules to combine in a single InvertedIndex.
			 */
			unsigned short blocksize_;
			
			
			/**
			 * Size of common counts matrix. Depends on blocksize_.
			 */
			LongSize cc_matrix_size_;
			
			
			/**
			 * Similarity cutoff.
			 */
			float cutoff_;
			
			
			/**
			 * Floating point precision.
			 */
			float precision_;
			
			
			/**
			 * If true, nearest neighbours are stored during connected components similarity calculations.
			 */
			bool store_nns_;
			
			
			/**
			 * Shared counter which initially stores the number of pairwise comparisons to be performed.
			 * This value is used to map into cells, i.e. rows and columns of a matrices.
			 * This value is decreased during calculations and finally 0 to indicate that all calculations have been performed!
			 */
			LongSize n_comparisons_;
			LongSize n_comparisons_backup_;
			
			
			/**
			 * Verbosity level.
			 */
			int verbosity_;
			
			
			/**
			 * Pointers to all leaf clusters.
			 */
			std::vector<Cluster*> leaf_clusters_;
			
			
			/**
			 * Map to store all internal clusters with same average linkage similarity.
			 */
			std::map<float, std::vector<Cluster*> > internal_clusters_;
			
			
			/**
			 * Active clusters. Used for various purposes.
			 */
			std::vector<Cluster*> vec_actives_;
			
			/**
			 * Inactive clusters. Used during RNN Parallel algorithm.
			 */
			std::vector<Cluster*> vec_inactives_;
			
			
			/**
			 * Similarity matrix stored as 1D array. Used for various purposes.
			 */
			float* sim_matrix_;
			
			/**
			 * Double precision similarity matrix stored as 1D array. Used for various purposes.
			 */
			double* dprec_sim_matrix_;
			
			
			/**
			 * InvertedIndices of all active clusters. Used during RNN Parallel algorithm.
			 */
			std::vector<std::pair<std::vector<InvertedIndex*>, unsigned int> > active_iids_;
			
			
			/**
			 * Number of active clusters to be gouped together for InvertedIndex generation.
			 * This is used to prevent InvertedIndices of size 1 in the early stage of the RNN Parallel algorithm.
			 */
			unsigned int active_iids_size_;
			
			
			/**
			 * Foremost cluster in Nearest Neighbour Chain.
			 */
			Cluster* nn_chain_tip_;
			
			
			/**
			 * Current size of the Nearest Neighbour Chain.
			 */
			unsigned int nn_chain_size_;
			
			
			/**
			 * Iterator pointing to current nearest neighbour of nn_chain_tip_.
			 */
			std::vector<Cluster*>::iterator current_nn_;
			
			
			/**
			 * Similarity of nn_chain_tip to current nearest neighbour:
			 */
			float current_nn_sim_;
			
			
			/**
			 * Clustering method to use.
			 */
			unsigned short clustering_method_;
			
			
			/**
			 * Total number of already created clusters.
			 */
			unsigned int n_clusters_;
			
			
			/**
			 * Maximal number of clusters / molecules for which Nearest Neighbour Chain method can be applied.
			 * THIS VALUE DETERMINES THE SIZE OF THE SIMILARITY MATRIX WHICH GROWS QUADRATIC. THUS THIS VALUE HAS TO BE HANDLED WITH CARE.
			 */
			unsigned int max_clusters_;
			
			
			/**
			 * Setup routine.
			 * @param options User defined options to overwrite default settings.
			 */
			void setup(const Options& options);
			
			
			/**
			 * Cleaning up routine.
			 * Destroys allocated memory.
			 */
			void clear();
			
			
			/**
			 * Assignment routine.
			 * @param bfm Class to be used for assignment.
			 */
			void assign(const BinaryFingerprintMethods& bfm);
			
			
			/**
			 * Check if input data (lib_features_ and selection) are valid. 
			 * @param selection Indices to molecule feature vectors in lib_features_ which should be checked.
			 * @return Returns true if data is consistent, false if there are inconsistencies.
			 */
			bool checkInputData(const std::vector<unsigned int>& selection) const;
			
			
			/**
			 * Allocates memory for thread data structures.
			 * @param blocksize Blocksize of InvertedIndices.
			 * @param dataset_size Determines the size of integer and floating point arrays for various purposes.
			 * @param active_iids_size Determines the size of floating or double precision matrices which are used to store shared feature counts or sum up similarities.
			 */
			void createThreadData(const unsigned int blocksize, const unsigned int dataset_size, const unsigned int active_iids_size);
			
			
			/**
			 * Deallocates memory from thread data structs.
			 */
			void destroyThreadData();
			
			
			/**
			 * Creates a single InvertedIndex from a set of molecules.
			 * @param members Pairs of pointers to feature vectors in lib_features_ or query_features_ and the ID of the cluster to which the molecule belongs.
			 * @return Pointer to the created InvertedIndex.
			 */
			InvertedIndex* createInvertedIndex(const std::vector<std::pair<const std::vector<unsigned short>*, unsigned int> >& members);
			
			
			/**
			 * Destroys an InvertedIndex.
			 * @param ii InvertedIndex to be destroyed.
			 */
			void destroyInvertedIndex(InvertedIndex* ii);
			
			
			/**
			 * Destroy all InvertedIndices in a vector of InvertedIndices.
			 * @param ii_destroy Vector with InvertedIndices to be destroyed. Vector is finally empty.
			 */
			void destroyInvertedIndices(InvertedIndices& ii_destroy);
			
			
			/**
			 * Create InvertedIndices of same size for a set of input molecules.
			 * @param molecules Feature vectors of molecules for which the InvertedIndices are created.
			 * @param target Reference to vector of InvertedIndices to store the created InvertedIndices (either lib_blocks or query_blocks).
			 */
			void createInvertedIndices(const std::vector<std::pair<const std::vector<unsigned short>*, unsigned int> >& molecules, InvertedIndices& ii_target);
			
			
			/**
			 * Set blocksize. As a side effect, the size of cc_matrix_size_ is set appropriately.
			 * @param blocksize Blocksize to be set.
			 */
			void setBlockSize(const unsigned short blocksize);
			
			
			/**
			 * Calculates row and column indices for an UPPER TRIANGULAR matrix (includes diagonal) from an index into an 1D array.
			 * @param row Reference to return the calculated row index of the matrix cell [0, n-1].
			 * @param col Reference to return the calculated column index of the matrix cell [0, n-1].
			 * @param array_index The index of the cell in the 1D array [0, ((n*n+n)/2)-1].
			 */
			void arrayToUpperTriangluarMatrix(LongSize& row, LongSize& col, const LongSize array_index) const;
			
			
			/**
			 * Calculates row and column indices for a STRICT UPPER TRIANGULAR matrix (excludes diagonal) from an index into an 1D array.
			 * @param row Reference to return the calculated row index of the matrix cell [0, n-1].
			 * @param col Reference to return the calculated column index of the matrix cell [0, n-1].
			 * @param array_index The index of the cell in the 1D array [0, ((n*n-n)/2)-1].
			 */
			void arrayToStrictUpperTriangularMatrix(LongSize& row, LongSize& col, const LongSize array_index) const;
			 
			
			/**
			 * Calculates 1D array index from row and column indices of a STRICT UPPER TRIANGULAR matrix (excludes diagonal).
			 * @param row Row index of the matrix cell [0, n-1].
			 * @param col Column index of the matrix cell [0, n-1].
			 * @return Array_index index of the cell in the 1D array [0, ((n*n-n)/2)-1].
			 */
			LongSize strictUpperTriangularMatrixToArray(const LongSize row, const LongSize col) const;
			
			
			/**
			 * Routine to coordiante a shared counter n_comparisons_ which is used to distribute the number of pairwise comparisons over parallel threads.
			 * @param index Integer reference to return index value.
			 * @return True if index is valid. False if n_comparisons_ = 0. Thus no comparisons are left.
			 */
			bool getNextComparisonIndex(LongSize& index);
			
			
			/**
			 * Compare similarity and ID of two items (molecules or clusters).
			 * @param a_sim Similarity of element a to be compared.
			 * @param b_sim Similarity of element b to be compared.
			 * @param a_id ID of element a to be compared.
			 * @param b_id ID of element b to be compared.
			 * @return True if: (a_sim > b_sim) || (a_sim == b_sim && a_id < b_id). Otherwise false.
			 */
			bool checkSimilaritySwitch(const float a_sim, const float b_sim, const unsigned int a_id, const unsigned int b_id) const;
			
			
			/**
			 * Shared feature count calculation: InvertedIndices ii1 and ii2 are both size 1.
			 * @param ii1 FeatureList of first InvertedIndex to be compared.
			 * @param ii2 FeatureList of second InvertedIndex to be compared.
			 * @param cc_count Counter which finally stores the number of shared features.
			 */
			void calculateCommonCounts_1_1(const FeatureList* ii1, const FeatureList* ii2, unsigned short& cc_count);
			
			
			/**
			 * Shared feature count calculation: InvertedIndices ii1 has size 1.
			 * @param ii1 FeatureList of first InvertedIndex to be compared.
			 * @param ii2 FeatureList of second InvertedIndex to be compared.
			 * @param cc_row 1D integer array which finally stores the number of shared features of (ii_1 X ii_2).
			 */
			void calculateCommonCounts_1_N(const FeatureList* ii1, const FeatureList* ii2, unsigned short* cc_row);
			
			
			/**
			 * Shared feature count calculation: InvertedIndices ii1 and ii2 have variable size.
			 * @param ii1 FeatureList of first InvertedIndex to be compared.
			 * @param ii2 FeatureList of second InvertedIndex to be compared.
			 * @param cc_matrix 2D integer array which finally stores the number of shared features of (ii_1 X ii_2).
			 */
			void calculateCommonCounts_M_N(const InvertedIndex* ii_1, const InvertedIndex* ii_2, unsigned short** cc_matrix);
			
			
			/**
			 * Calculation of similarity coefficients and writing of similarities above cutoff to outfile.
			 * @param query_index Position of InvertedIndex in query_iindices_.
			 * @param lib_index Position of InvertedIndex in lib_iindices_.
			 * @param cc_matrix Shared feature counts matrix between query_index and lib_index.
			 * @param outfile Reference to output file stream to write result values.
			 */
			void cutoffSearchSimilarities(const unsigned int query_index, const unsigned int lib_index, unsigned short** cc_matrix, File& outfile);
			
			
			/**
			 * Thread routine to run similarity calculation.
			 * @param thread_id ID if corresponding thread.
			 */
			void cutoffSearchThread(const unsigned int thread_id);
			
			
			/**
			 * Function pointer to access the different functions for pairwise similarity calculation.
			 * @param ii1_index Position of first InvertedIndex to be compared.
			 * @param ii2_index Position of second InvertedIndex to be compared.
			 * @param t_data Pointer to thread data structure.
			 */
			typedef void (BinaryFingerprintMethods::*PairwiseSimilaritiesBase)(const unsigned int ii1_index, const unsigned int ii2_index, ThreadData* t_data);
			PairwiseSimilaritiesBase pairwiseSimilaritiesBase;
			
			
			/**
			 * Pairwise similarity calculation which implements nearest neighbour search.
			 * @param ii1_index Position of first InvertedIndex to be compared.
			 * @param ii2_index Position of second InvertedIndex to be compared.
			 * @param t_data Pointer to thread data structure.
			 */
			void pairwiseSimilaritiesNearestNeighbours(const unsigned int ii1_index, const unsigned int ii2_index, ThreadData* t_data);
			
			
			/**
			 * Pairwise similarity calculation which implements calculation of a similarity matrix.
			 * @param ii1_index Position of first InvertedIndex to be compared.
			 * @param ii2_index Position of second InvertedIndex to be compared.
			 * @param t_data Pointer to thread data structure.
			 */
			void pairwiseSimilaritiesStoredMatrix(const unsigned int ii1_index, const unsigned int ii2_index, ThreadData* t_data);
			
			
			/**
			 * Pairwise similarity calculation which implements similarity network generation for connected components decomposition.
			 * If store_nns_ is true, nearest neighbours are also calculated.
			 * @param ii1_index Position of first InvertedIndex to be compared.
			 * @param ii2_index Position of second InvertedIndex to be compared.
			 * @param t_data Pointer to thread data structure.
			 */
			void pairwiseSimilaritiesConnectedComponents(const unsigned int ii1_index, const unsigned int ii2_index, ThreadData* t_data);
			
			
			/**
			 * Pairwise similarity calculation which implements summing up the similarity of every compound to all others. Used for medoid calculation.
			 * @param ii1_index Position of first InvertedIndex to be compared.
			 * @param ii2_index Position of second InvertedIndex to be compared.
			 * @param t_data Pointer to thread data structure.
			 */
			void pairwiseSimilaritiesMedoids(const unsigned int ii1_index, const unsigned int ii2_index, ThreadData* t_data);
			
			
			/**
			 * Base function which runs all pairwise comparisons of molecules indexed by selection.
			 * @param selection Indices of feature vectors in lib_features_ for which pairwise comparisons should be performed.
			 * @param nn_data If store_nns_ is true, nn_data will finally store nearest neighbour information for every molecule.
			 * @return True if calculation was successful.
			 */
			bool pairwiseSimilarities(const std::vector<unsigned int>& selection, std::vector<std::pair<unsigned int, float> >& nn_data);
			
			
			/**
			 * Thread function for all pairwise comparisons.
			 * @param thread_id ID of corresponding BOOST thread.
			 */
			void pairwiseSimilaritiesThread(const unsigned int thread_id);
			
			
			/**
			 * Similarity calculations of active clusters. Used by RNN parallel algorithm.
			 * @param ii1 Pointer to first InvertedIndex to be compared.
			 * @param ii1 Pointer to second InvertedIndex to be compared.
			 * @param cc_matrix Integer matrix which stores shared features counts for all pairs of molecules in (ii1 X ii2).
			 * @param sim_matrix Double precision matrix to sum up iter-cluster similarities of two disjunct clusters.
			 */
			void calculateParallelSimilaritiesActives(const InvertedIndex* ii1, const InvertedIndex* ii2, unsigned short** cc_matrix, double** sim_matrix);
			
			
			/**
			 * Thread function for pairwise similarity calculations of active clusters. Used by RNN parallel algorithm.
			 * @param thread_id ID of corresponding BOOST thread.
			 */
			void calculateParallelSimilaritiesActivesThread(const unsigned int thread_id);
			
			
			/**
			 * Similarity calculations of active clusters against inactive clusters. Used by RNN parallel algorithm.
			 * @param ii1 Pointer to first InvertedIndex to be compared.
			 * @param ii1 Pointer to second InvertedIndex to be compared.
			 * @param cc_matrix Integer matrix which stores shared features counts for all pairs of molecules in (ii1 X ii2).
			 * @param sim_matrix Double precision matrix to sum up iter-cluster similarities of two disjunct clusters.
			 */
			void calculateParallelSimilarities(const InvertedIndex* ii1, const InvertedIndex* ii2, unsigned short** cc_matrix, double** sim_matrix);
			
			
			/**
			 * Thread function for pairwise similarity calculations of active clusters aginst inactive clusters. Used by RNN parallel algorithm.
			 * @param thread_id ID of corresponding BOOST thread.
			 */
			void calculateParallelSimilaritiesThread(const unsigned int thread_id);
			
			
			/**
			 * Calculate similarity of newly merged cluster to another cluster accodrding to Lance-Williams update formula.
			 * G. Lance and W. Williams, Comput J. (1967), 9, 373-380. (doi: 10.1093/comjnl/9.4.373).
			 * @param merged_cluster Pointer to newly merged cluster
			 * @param current Pointer to cluster to which similarity is calculated.
			 */
			void similarityUpdateAverageLinkage(const Cluster* merged_cluster, const Cluster* current);
			
			/**
			 * Thread function to update similarities of newly merged cluster to another cluster accodrding to Lance-Williams update formula.
			 * @param thread_id ID of corresponding BOOST thread.
			 * @param merged_cluster Pointer to newly merged cluster.
			 */
			void similarityUpdateAverageLinkageThread(const unsigned int thread_id, const Cluster* merged_cluster);
			
			
			/**
			 * Summing up the inter-cluster similarities of two InvertedIndices where ii1 has size 1.
			 * @param ii1 Pointer to first InvertedIndex to be compared which has size 1. 
			 * @param ii2 Pointer to second InvertedIndex to be compared.
			 * @param cc_row Integer array which stores shared feature counts.
			 * @param sim_sum Variable to sum up inter-cluster similarities.
			 */
			void clusterSimilaritySum_1_N(const InvertedIndex* ii1, const InvertedIndex* ii2, const unsigned short* cc_row, double& sim_sum);
			
			
			/**
			 * Summing up the inter-cluster similarities of two InvertedIndices.
			 * @param ii1 Pointer to first InvertedIndex to be compared.
			 * @param ii2 Pointer to second InvertedIndex to be compared.
			 * @param cc_matrix Integer matrix which stores shared feature counts.
			 * @param sim_sum Variable to sum up inter-cluster similarities.
			 */
			void clusterSimilaritySum_M_N(const InvertedIndex* ii1, const InvertedIndex* ii2, unsigned short** cc_matrix, double& sim_sum);
			
			
			/**
			 * Calculates similarity matrix from all active clusters. Used during switch from RNN parallel to Nearest Neighbour Chain algorithm.
			 * @param thread_id ID of corresponding BOOST thread.
			 */
			void similarityMatrixFromClustersThread(const unsigned int thread_id);
			
			
			/**
			 * Clustering algorithm using the parallel RNN approach described by F. Murtagh.\n
			 * F. Murtagh, Compstat Lectures vol. 4, 1985, Physica-Verlag Würzburg-Wien. (Volume titel: Multidimensional clustering algorithms)
			 * @param root Finally points to the root node, i.e. the final cluster comprising all other nodes.
			 */
			void averageLinkageParallel(Cluster*& root);
			
			
			/**
			 * Nearest Neighbour Chain (NNChain) algorithm initially described by F. Murtagh.\n
			 * F. Murtagh, Compstat Lectures vol. 4, 1985, Physica-Verlag Würzburg-Wien. (Volume titel: Multidimensional clustering algorithms)
			 * @param root Finally points to the root node, i.e. the final cluster comprising all other nodes.
			 */
			void NNChainCore(Cluster*& root);
			
			
			/**
			 * Method calculates the Kelley-Gardner-Sutcliffe penalties for cluster level selection.\n
			 * L.A. Kelley, S.P. Gardner and M.J. Sutcliffe, Protein Eng. (1996), 9, 1063-1065. (doi: 10.1093/protein/9.11.1063).
			 * @param cluster_selection Vector to store the final cluster selection.
			 */
			void clusterSelectionKGS(std::map<unsigned int, std::vector<unsigned int> >& cluster_selection);
			
			
			/**
			 * Recursively walk down the dendrogram and assign cluster IDs to leaf nodes (i.e. molecules).
			 * @param cl Pointer to urrent cluster.
			 * @param cluster_id Cluster ID of this subtree.
			 */
			void enumerateClusterMembers(Cluster* cl, unsigned int cluster_id);
			
			
			/**
			 * Nearest Neighbour Chain Algorithm: Calculate the nearest neighbour of the nn_chain_tip_ within the set of active clusters.
			 */
			void nextNearestNeighbour();
			
			
			/**
			 * Nearest Neighbour Chain Algorithm: Append nearest neighbour of current the nn_chain_tip_ to Nearest Neighbour Chain.
			 */
			void moveNearestNeighbour();
			
			
			/**
			 * Merge two clusters.
			 * @param c1 Pointer to first cluster to merge.
			 * @param c2 Pointer to second cluster to merge.
			 * @param sim_sum Similarity sum of inter-cluster similarities of c1 and c2 (Nearest Neighbour Chain algorithm) or similarity of c1 to c2 (RNN parallel).
			 * @return Pointer to merged and newly created cluster.
			 */
			Cluster* mergeClusters(Cluster* c1, Cluster* c2, double sim_sum);
			
			
			/**
			 * Function creates a new cluster.
			 * @return Pointer to newly formed cluster.
			 */
			Cluster* createCluster();
			
			/**
			 * Switch from RNN parallel stored data to Nearest Neighbour Chain stored matrix algorithm.
			 */
			void switchStorageMethod();
			
			
			/**
			 * Clean up when finished clustering.
			 */
			void finalizeClustering();
	};
} // namespace BALL 

#endif // BALL_STRUCTURE_BINARYFINGERPRINTMETHODS_H