//
// C++ standalone verion of fastcluster by Daniel Müllner
//
// Copyright: Christoph Dalitz, 2018
// Daniel Müllner, 2011
// License: BSD style license
// (see the file LICENSE for details)
//
# include <vector>
# include <algorithm>
# include <cmath>
extern " C " {
# include "fastcluster.h"
}
// Code by Daniel Müllner
// workaround to make it usable as a standalone version (without R)
bool fc_isnan ( double x ) { return false ; }
# include "fastcluster_dm.cpp"
# include "fastcluster_R_dm.cpp"
extern " C " {
//
// Assigns cluster labels (0, ..., nclust-1) to the n points such
// that the cluster result is split into nclust clusters.
//
// Input arguments:
// n = number of observables
// merge = clustering result in R format
// nclust = number of clusters
// Output arguments:
// labels = allocated integer array of size n for result
//
void cutree_k ( int n , const int * merge , int nclust , int * labels ) {
int k , m1 , m2 , j , l ;
if ( nclust > n | | nclust < 2 ) {
for ( j = 0 ; j < n ; j + + ) labels [ j ] = 0 ;
return ;
}
// assign to each observable the number of its last merge step
// beware: indices of observables in merge start at 1 (R convention)
std : : vector < int > last_merge ( n , 0 ) ;
for ( k = 1 ; k < = ( n - nclust ) ; k + + ) {
// (m1,m2) = merge[k,]
m1 = merge [ k - 1 ] ;
m2 = merge [ n - 1 + k - 1 ] ;
if ( m1 < 0 & & m2 < 0 ) { // both single observables
last_merge [ - m1 - 1 ] = last_merge [ - m2 - 1 ] = k ;
}
else if ( m1 < 0 | | m2 < 0 ) { // one is a cluster
if ( m1 < 0 ) { j = - m1 ; m1 = m2 ; } else j = - m2 ;
// merging single observable and cluster
for ( l = 0 ; l < n ; l + + )
if ( last_merge [ l ] = = m1 )
last_merge [ l ] = k ;
last_merge [ j - 1 ] = k ;
}
else { // both cluster
for ( l = 0 ; l < n ; l + + ) {
if ( last_merge [ l ] = = m1 | | last_merge [ l ] = = m2 )
last_merge [ l ] = k ;
}
}
}
// assign cluster labels
int label = 0 ;
std : : vector < int > z ( n , - 1 ) ;
for ( j = 0 ; j < n ; j + + ) {
if ( last_merge [ j ] = = 0 ) { // still singleton
labels [ j ] = label + + ;
} else {
if ( z [ last_merge [ j ] ] < 0 ) {
z [ last_merge [ j ] ] = label + + ;
}
labels [ j ] = z [ last_merge [ j ] ] ;
}
}
}
//
// Assigns cluster labels (0, ..., nclust-1) to the n points such
// that the hierarchical clustering is stopped when cluster distance >= cdist
//
// Input arguments:
// n = number of observables
// merge = clustering result in R format
// height = cluster distance at each merge step
// cdist = cutoff cluster distance
// Output arguments:
// labels = allocated integer array of size n for result
//
void cutree_cdist ( int n , const int * merge , double * height , double cdist , int * labels ) {
int k ;
for ( k = 0 ; k < ( n - 1 ) ; k + + ) {
if ( height [ k ] > = cdist ) {
break ;
}
}
cutree_k ( n , merge , n - k , labels ) ;
}
//
// Hierarchical clustering with one of Daniel Muellner's fast algorithms
//
// Input arguments:
// n = number of observables
// distmat = condensed distance matrix, i.e. an n*(n-1)/2 array representing
// the upper triangle (without diagonal elements) of the distance
// matrix, e.g. for n=4:
// d00 d01 d02 d03
// d10 d11 d12 d13 -> d01 d02 d03 d12 d13 d23
// d20 d21 d22 d23
// d30 d31 d32 d33
// method = cluster metric (see enum method_code)
// Output arguments:
// merge = allocated (n-1)x2 matrix (2*(n-1) array) for storing result.
// Result follows R hclust convention:
// - observabe indices start with one
// - merge[i][] contains the merged nodes in step i
// - merge[i][j] is negative when the node is an atom
// height = allocated (n-1) array with distances at each merge step
// Return code:
// 0 = ok
// 1 = invalid method
//
int hclust_fast ( int n , double * distmat , int method , int * merge , double * height ) {
// call appropriate culstering function
cluster_result Z2 ( n - 1 ) ;
if ( method = = HCLUST_METHOD_SINGLE ) {
// single link
MST_linkage_core ( n , distmat , Z2 ) ;
}
else if ( method = = HCLUST_METHOD_COMPLETE ) {
// complete link
NN_chain_core < METHOD_METR_COMPLETE , t_float > ( n , distmat , NULL , Z2 ) ;
}
else if ( method = = HCLUST_METHOD_AVERAGE ) {
// best average distance
double * members = new double [ n ] ;
for ( int i = 0 ; i < n ; i + + ) members [ i ] = 1 ;
NN_chain_core < METHOD_METR_AVERAGE , t_float > ( n , distmat , members , Z2 ) ;
delete [ ] members ;
}
else if ( method = = HCLUST_METHOD_MEDIAN ) {
// best median distance (beware: O(n^3))
generic_linkage < METHOD_METR_MEDIAN , t_float > ( n , distmat , NULL , Z2 ) ;
}
else if ( method = = HCLUST_METHOD_CENTROID ) {
// best centroid distance (beware: O(n^3))
double * members = new double [ n ] ;
for ( int i = 0 ; i < n ; i + + ) members [ i ] = 1 ;
generic_linkage < METHOD_METR_CENTROID , t_float > ( n , distmat , members , Z2 ) ;
delete [ ] members ;
}
else {
return 1 ;
}
int * order = new int [ n ] ;
if ( method = = HCLUST_METHOD_MEDIAN | | method = = HCLUST_METHOD_CENTROID ) {
generate_R_dendrogram < true > ( merge , height , order , Z2 , n ) ;
} else {
generate_R_dendrogram < false > ( merge , height , order , Z2 , n ) ;
}
delete [ ] order ; // only needed for visualization
return 0 ;
}
// Build condensed distance matrix
// Input arguments:
// n = number of observables
// m = dimension of observable
// Output arguments:
// out = allocated integer array of size n * (n - 1) / 2 for result
void hclust_pdist ( int n , int m , double * pts , double * out ) {
int ii = 0 ;
for ( int i = 0 ; i < n ; i + + ) {
for ( int j = i + 1 ; j < n ; j + + ) {
// Compute euclidian distance
double d = 0 ;
for ( int k = 0 ; k < m ; k + + ) {
double error = pts [ i * m + k ] - pts [ j * m + k ] ;
d + = ( error * error ) ;
}
out [ ii ] = d ; //sqrt(d);
ii + + ;
}
}
}
void cluster_points_centroid ( int n , int m , double * pts , double dist , int * idx ) {
double * pdist = new double [ n * ( n - 1 ) / 2 ] ;
int * merge = new int [ 2 * ( n - 1 ) ] ;
double * height = new double [ n - 1 ] ;
hclust_pdist ( n , m , pts , pdist ) ;
hclust_fast ( n , pdist , HCLUST_METHOD_CENTROID , merge , height ) ;
cutree_cdist ( n , merge , height , dist , idx ) ;
delete [ ] pdist ;
delete [ ] merge ;
delete [ ] height ;
}
}
