You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1794 lines
51 KiB
1794 lines
51 KiB
/*
|
|
fastcluster: Fast hierarchical clustering routines for R and Python
|
|
|
|
Copyright © 2011 Daniel Müllner
|
|
<http://danifold.net>
|
|
|
|
This library implements various fast algorithms for hierarchical,
|
|
agglomerative clustering methods:
|
|
|
|
(1) Algorithms for the "stored matrix approach": the input is the array of
|
|
pairwise dissimilarities.
|
|
|
|
MST_linkage_core: single linkage clustering with the "minimum spanning
|
|
tree algorithm (Rohlfs)
|
|
|
|
NN_chain_core: nearest-neighbor-chain algorithm, suitable for single,
|
|
complete, average, weighted and Ward linkage (Murtagh)
|
|
|
|
generic_linkage: generic algorithm, suitable for all distance update
|
|
formulas (Müllner)
|
|
|
|
(2) Algorithms for the "stored data approach": the input are points in a
|
|
vector space.
|
|
|
|
MST_linkage_core_vector: single linkage clustering for vector data
|
|
|
|
generic_linkage_vector: generic algorithm for vector data, suitable for
|
|
the Ward, centroid and median methods.
|
|
|
|
generic_linkage_vector_alternative: alternative scheme for updating the
|
|
nearest neighbors. This method seems faster than "generic_linkage_vector"
|
|
for the centroid and median methods but slower for the Ward method.
|
|
|
|
All these implementation treat infinity values correctly. They throw an
|
|
exception if a NaN distance value occurs.
|
|
*/
|
|
|
|
// Older versions of Microsoft Visual Studio do not have the fenv header.
|
|
#ifdef _MSC_VER
|
|
#if (_MSC_VER == 1500 || _MSC_VER == 1600)
|
|
#define NO_INCLUDE_FENV
|
|
#endif
|
|
#endif
|
|
// NaN detection via fenv might not work on systems with software
|
|
// floating-point emulation (bug report for Debian armel).
|
|
#ifdef __SOFTFP__
|
|
#define NO_INCLUDE_FENV
|
|
#endif
|
|
#ifdef NO_INCLUDE_FENV
|
|
#pragma message("Do not use fenv header.")
|
|
#else
|
|
#pragma message("Use fenv header. If there is a warning about unknown #pragma STDC FENV_ACCESS, this can be ignored.")
|
|
#pragma STDC FENV_ACCESS on
|
|
#include <fenv.h>
|
|
#endif
|
|
|
|
#include <cmath> // for std::pow, std::sqrt
|
|
#include <cstddef> // for std::ptrdiff_t
|
|
#include <limits> // for std::numeric_limits<...>::infinity()
|
|
#include <algorithm> // for std::fill_n
|
|
#include <stdexcept> // for std::runtime_error
|
|
#include <string> // for std::string
|
|
|
|
#include <cfloat> // also for DBL_MAX, DBL_MIN
|
|
#ifndef DBL_MANT_DIG
|
|
#error The constant DBL_MANT_DIG could not be defined.
|
|
#endif
|
|
#define T_FLOAT_MANT_DIG DBL_MANT_DIG
|
|
|
|
#ifndef LONG_MAX
|
|
#include <climits>
|
|
#endif
|
|
#ifndef LONG_MAX
|
|
#error The constant LONG_MAX could not be defined.
|
|
#endif
|
|
#ifndef INT_MAX
|
|
#error The constant INT_MAX could not be defined.
|
|
#endif
|
|
|
|
#ifndef INT32_MAX
|
|
#ifdef _MSC_VER
|
|
#if _MSC_VER >= 1600
|
|
#define __STDC_LIMIT_MACROS
|
|
#include <stdint.h>
|
|
#else
|
|
typedef __int32 int_fast32_t;
|
|
typedef __int64 int64_t;
|
|
#endif
|
|
#else
|
|
#define __STDC_LIMIT_MACROS
|
|
#include <stdint.h>
|
|
#endif
|
|
#endif
|
|
|
|
#define FILL_N std::fill_n
|
|
#ifdef _MSC_VER
|
|
#if _MSC_VER < 1600
|
|
#undef FILL_N
|
|
#define FILL_N stdext::unchecked_fill_n
|
|
#endif
|
|
#endif
|
|
|
|
// Suppress warnings about (potentially) uninitialized variables.
|
|
#ifdef _MSC_VER
|
|
#pragma warning (disable:4700)
|
|
#endif
|
|
|
|
#ifndef HAVE_DIAGNOSTIC
|
|
#if __GNUC__ > 4 || (__GNUC__ == 4 && (__GNUC_MINOR__ >= 6))
|
|
#define HAVE_DIAGNOSTIC 1
|
|
#endif
|
|
#endif
|
|
|
|
#ifndef HAVE_VISIBILITY
|
|
#if __GNUC__ >= 4
|
|
#define HAVE_VISIBILITY 1
|
|
#endif
|
|
#endif
|
|
|
|
/* Since the public interface is given by the Python respectively R interface,
|
|
* we do not want other symbols than the interface initalization routines to be
|
|
* visible in the shared object file. The "visibility" switch is a GCC concept.
|
|
* Hiding symbols keeps the relocation table small and decreases startup time.
|
|
* See http://gcc.gnu.org/wiki/Visibility
|
|
*/
|
|
#if HAVE_VISIBILITY
|
|
#pragma GCC visibility push(hidden)
|
|
#endif
|
|
|
|
typedef int_fast32_t t_index;
|
|
#ifndef INT32_MAX
|
|
#define MAX_INDEX 0x7fffffffL
|
|
#else
|
|
#define MAX_INDEX INT32_MAX
|
|
#endif
|
|
#if (LONG_MAX < MAX_INDEX)
|
|
#error The integer format "t_index" must not have a greater range than "long int".
|
|
#endif
|
|
#if (INT_MAX > MAX_INDEX)
|
|
#error The integer format "int" must not have a greater range than "t_index".
|
|
#endif
|
|
typedef double t_float;
|
|
|
|
/* Method codes.
|
|
|
|
These codes must agree with the METHODS array in fastcluster.R and the
|
|
dictionary mthidx in fastcluster.py.
|
|
*/
|
|
enum method_codes {
|
|
// non-Euclidean methods
|
|
METHOD_METR_SINGLE = 0,
|
|
METHOD_METR_COMPLETE = 1,
|
|
METHOD_METR_AVERAGE = 2,
|
|
METHOD_METR_WEIGHTED = 3,
|
|
METHOD_METR_WARD = 4,
|
|
METHOD_METR_WARD_D = METHOD_METR_WARD,
|
|
METHOD_METR_CENTROID = 5,
|
|
METHOD_METR_MEDIAN = 6,
|
|
METHOD_METR_WARD_D2 = 7,
|
|
|
|
MIN_METHOD_CODE = 0,
|
|
MAX_METHOD_CODE = 7
|
|
};
|
|
|
|
enum method_codes_vector {
|
|
// Euclidean methods
|
|
METHOD_VECTOR_SINGLE = 0,
|
|
METHOD_VECTOR_WARD = 1,
|
|
METHOD_VECTOR_CENTROID = 2,
|
|
METHOD_VECTOR_MEDIAN = 3,
|
|
|
|
MIN_METHOD_VECTOR_CODE = 0,
|
|
MAX_METHOD_VECTOR_CODE = 3
|
|
};
|
|
|
|
// self-destructing array pointer
|
|
template <typename type>
|
|
class auto_array_ptr{
|
|
private:
|
|
type * ptr;
|
|
auto_array_ptr(auto_array_ptr const &); // non construction-copyable
|
|
auto_array_ptr& operator=(auto_array_ptr const &); // non copyable
|
|
public:
|
|
auto_array_ptr()
|
|
: ptr(NULL)
|
|
{ }
|
|
template <typename index>
|
|
auto_array_ptr(index const size)
|
|
: ptr(new type[size])
|
|
{ }
|
|
template <typename index, typename value>
|
|
auto_array_ptr(index const size, value const val)
|
|
: ptr(new type[size])
|
|
{
|
|
FILL_N(ptr, size, val);
|
|
}
|
|
~auto_array_ptr() {
|
|
delete [] ptr; }
|
|
void free() {
|
|
delete [] ptr;
|
|
ptr = NULL;
|
|
}
|
|
template <typename index>
|
|
void init(index const size) {
|
|
ptr = new type [size];
|
|
}
|
|
template <typename index, typename value>
|
|
void init(index const size, value const val) {
|
|
init(size);
|
|
FILL_N(ptr, size, val);
|
|
}
|
|
inline operator type *() const { return ptr; }
|
|
};
|
|
|
|
struct node {
|
|
t_index node1, node2;
|
|
t_float dist;
|
|
};
|
|
|
|
inline bool operator< (const node a, const node b) {
|
|
return (a.dist < b.dist);
|
|
}
|
|
|
|
class cluster_result {
|
|
private:
|
|
auto_array_ptr<node> Z;
|
|
t_index pos;
|
|
|
|
public:
|
|
cluster_result(const t_index size)
|
|
: Z(size)
|
|
, pos(0)
|
|
{}
|
|
|
|
void append(const t_index node1, const t_index node2, const t_float dist) {
|
|
Z[pos].node1 = node1;
|
|
Z[pos].node2 = node2;
|
|
Z[pos].dist = dist;
|
|
++pos;
|
|
}
|
|
|
|
node * operator[] (const t_index idx) const { return Z + idx; }
|
|
|
|
/* Define several methods to postprocess the distances. All these functions
|
|
are monotone, so they do not change the sorted order of distances. */
|
|
|
|
void sqrt() const {
|
|
for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
|
|
ZZ->dist = std::sqrt(ZZ->dist);
|
|
}
|
|
}
|
|
|
|
void sqrt(const t_float) const { // ignore the argument
|
|
sqrt();
|
|
}
|
|
|
|
void sqrtdouble(const t_float) const { // ignore the argument
|
|
for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
|
|
ZZ->dist = std::sqrt(2*ZZ->dist);
|
|
}
|
|
}
|
|
|
|
#ifdef R_pow
|
|
#define my_pow R_pow
|
|
#else
|
|
#define my_pow std::pow
|
|
#endif
|
|
|
|
void power(const t_float p) const {
|
|
t_float const q = 1/p;
|
|
for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
|
|
ZZ->dist = my_pow(ZZ->dist,q);
|
|
}
|
|
}
|
|
|
|
void plusone(const t_float) const { // ignore the argument
|
|
for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
|
|
ZZ->dist += 1;
|
|
}
|
|
}
|
|
|
|
void divide(const t_float denom) const {
|
|
for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
|
|
ZZ->dist /= denom;
|
|
}
|
|
}
|
|
};
|
|
|
|
class doubly_linked_list {
|
|
/*
|
|
Class for a doubly linked list. Initially, the list is the integer range
|
|
[0, size]. We provide a forward iterator and a method to delete an index
|
|
from the list.
|
|
|
|
Typical use: for (i=L.start; L<size; i=L.succ[I])
|
|
or
|
|
for (i=somevalue; L<size; i=L.succ[I])
|
|
*/
|
|
public:
|
|
t_index start;
|
|
auto_array_ptr<t_index> succ;
|
|
|
|
private:
|
|
auto_array_ptr<t_index> pred;
|
|
// Not necessarily private, we just do not need it in this instance.
|
|
|
|
public:
|
|
doubly_linked_list(const t_index size)
|
|
// Initialize to the given size.
|
|
: start(0)
|
|
, succ(size+1)
|
|
, pred(size+1)
|
|
{
|
|
for (t_index i=0; i<size; ++i) {
|
|
pred[i+1] = i;
|
|
succ[i] = i+1;
|
|
}
|
|
// pred[0] is never accessed!
|
|
//succ[size] is never accessed!
|
|
}
|
|
|
|
~doubly_linked_list() {}
|
|
|
|
void remove(const t_index idx) {
|
|
// Remove an index from the list.
|
|
if (idx==start) {
|
|
start = succ[idx];
|
|
}
|
|
else {
|
|
succ[pred[idx]] = succ[idx];
|
|
pred[succ[idx]] = pred[idx];
|
|
}
|
|
succ[idx] = 0; // Mark as inactive
|
|
}
|
|
|
|
bool is_inactive(t_index idx) const {
|
|
return (succ[idx]==0);
|
|
}
|
|
};
|
|
|
|
// Indexing functions
|
|
// D is the upper triangular part of a symmetric (NxN)-matrix
|
|
// We require r_ < c_ !
|
|
#define D_(r_,c_) ( D[(static_cast<std::ptrdiff_t>(2*N-3-(r_))*(r_)>>1)+(c_)-1] )
|
|
// Z is an ((N-1)x4)-array
|
|
#define Z_(_r, _c) (Z[(_r)*4 + (_c)])
|
|
|
|
/*
|
|
Lookup function for a union-find data structure.
|
|
|
|
The function finds the root of idx by going iteratively through all
|
|
parent elements until a root is found. An element i is a root if
|
|
nodes[i] is zero. To make subsequent searches faster, the entry for
|
|
idx and all its parents is updated with the root element.
|
|
*/
|
|
class union_find {
|
|
private:
|
|
auto_array_ptr<t_index> parent;
|
|
t_index nextparent;
|
|
|
|
public:
|
|
union_find(const t_index size)
|
|
: parent(size>0 ? 2*size-1 : 0, 0)
|
|
, nextparent(size)
|
|
{ }
|
|
|
|
t_index Find (t_index idx) const {
|
|
if (parent[idx] != 0 ) { // a → b
|
|
t_index p = idx;
|
|
idx = parent[idx];
|
|
if (parent[idx] != 0 ) { // a → b → c
|
|
do {
|
|
idx = parent[idx];
|
|
} while (parent[idx] != 0);
|
|
do {
|
|
t_index tmp = parent[p];
|
|
parent[p] = idx;
|
|
p = tmp;
|
|
} while (parent[p] != idx);
|
|
}
|
|
}
|
|
return idx;
|
|
}
|
|
|
|
void Union (const t_index node1, const t_index node2) {
|
|
parent[node1] = parent[node2] = nextparent++;
|
|
}
|
|
};
|
|
|
|
class nan_error{};
|
|
#ifdef FE_INVALID
|
|
class fenv_error{};
|
|
#endif
|
|
|
|
static void MST_linkage_core(const t_index N, const t_float * const D,
|
|
cluster_result & Z2) {
|
|
/*
|
|
N: integer, number of data points
|
|
D: condensed distance matrix N*(N-1)/2
|
|
Z2: output data structure
|
|
|
|
The basis of this algorithm is an algorithm by Rohlf:
|
|
|
|
F. James Rohlf, Hierarchical clustering using the minimum spanning tree,
|
|
The Computer Journal, vol. 16, 1973, p. 93–95.
|
|
*/
|
|
t_index i;
|
|
t_index idx2;
|
|
doubly_linked_list active_nodes(N);
|
|
auto_array_ptr<t_float> d(N);
|
|
|
|
t_index prev_node;
|
|
t_float min;
|
|
|
|
// first iteration
|
|
idx2 = 1;
|
|
min = std::numeric_limits<t_float>::infinity();
|
|
for (i=1; i<N; ++i) {
|
|
d[i] = D[i-1];
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (d[i] < min) {
|
|
min = d[i];
|
|
idx2 = i;
|
|
}
|
|
else if (fc_isnan(d[i]))
|
|
throw (nan_error());
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
}
|
|
Z2.append(0, idx2, min);
|
|
|
|
for (t_index j=1; j<N-1; ++j) {
|
|
prev_node = idx2;
|
|
active_nodes.remove(prev_node);
|
|
|
|
idx2 = active_nodes.succ[0];
|
|
min = d[idx2];
|
|
for (i=idx2; i<prev_node; i=active_nodes.succ[i]) {
|
|
t_float tmp = D_(i, prev_node);
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (tmp < d[i])
|
|
d[i] = tmp;
|
|
else if (fc_isnan(tmp))
|
|
throw (nan_error());
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
if (d[i] < min) {
|
|
min = d[i];
|
|
idx2 = i;
|
|
}
|
|
}
|
|
for (; i<N; i=active_nodes.succ[i]) {
|
|
t_float tmp = D_(prev_node, i);
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (d[i] > tmp)
|
|
d[i] = tmp;
|
|
else if (fc_isnan(tmp))
|
|
throw (nan_error());
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
if (d[i] < min) {
|
|
min = d[i];
|
|
idx2 = i;
|
|
}
|
|
}
|
|
Z2.append(prev_node, idx2, min);
|
|
}
|
|
}
|
|
|
|
/* Functions for the update of the dissimilarity array */
|
|
|
|
inline static void f_single( t_float * const b, const t_float a ) {
|
|
if (*b > a) *b = a;
|
|
}
|
|
inline static void f_complete( t_float * const b, const t_float a ) {
|
|
if (*b < a) *b = a;
|
|
}
|
|
inline static void f_average( t_float * const b, const t_float a, const t_float s, const t_float t) {
|
|
*b = s*a + t*(*b);
|
|
#ifndef FE_INVALID
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (fc_isnan(*b)) {
|
|
throw(nan_error());
|
|
}
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
#endif
|
|
}
|
|
inline static void f_weighted( t_float * const b, const t_float a) {
|
|
*b = (a+*b)*.5;
|
|
#ifndef FE_INVALID
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (fc_isnan(*b)) {
|
|
throw(nan_error());
|
|
}
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
#endif
|
|
}
|
|
inline static void f_ward( t_float * const b, const t_float a, const t_float c, const t_float s, const t_float t, const t_float v) {
|
|
*b = ( (v+s)*a - v*c + (v+t)*(*b) ) / (s+t+v);
|
|
//*b = a+(*b)-(t*a+s*(*b)+v*c)/(s+t+v);
|
|
#ifndef FE_INVALID
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (fc_isnan(*b)) {
|
|
throw(nan_error());
|
|
}
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
#endif
|
|
}
|
|
inline static void f_centroid( t_float * const b, const t_float a, const t_float stc, const t_float s, const t_float t) {
|
|
*b = s*a - stc + t*(*b);
|
|
#ifndef FE_INVALID
|
|
if (fc_isnan(*b)) {
|
|
throw(nan_error());
|
|
}
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
#endif
|
|
}
|
|
inline static void f_median( t_float * const b, const t_float a, const t_float c_4) {
|
|
*b = (a+(*b))*.5 - c_4;
|
|
#ifndef FE_INVALID
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (fc_isnan(*b)) {
|
|
throw(nan_error());
|
|
}
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
#endif
|
|
}
|
|
|
|
template <method_codes method, typename t_members>
|
|
static void NN_chain_core(const t_index N, t_float * const D, t_members * const members, cluster_result & Z2) {
|
|
/*
|
|
N: integer
|
|
D: condensed distance matrix N*(N-1)/2
|
|
Z2: output data structure
|
|
|
|
This is the NN-chain algorithm, described on page 86 in the following book:
|
|
|
|
Fionn Murtagh, Multidimensional Clustering Algorithms,
|
|
Vienna, Würzburg: Physica-Verlag, 1985.
|
|
*/
|
|
t_index i;
|
|
|
|
auto_array_ptr<t_index> NN_chain(N);
|
|
t_index NN_chain_tip = 0;
|
|
|
|
t_index idx1, idx2;
|
|
|
|
t_float size1, size2;
|
|
doubly_linked_list active_nodes(N);
|
|
|
|
t_float min;
|
|
|
|
for (t_float const * DD=D; DD!=D+(static_cast<std::ptrdiff_t>(N)*(N-1)>>1);
|
|
++DD) {
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (fc_isnan(*DD)) {
|
|
throw(nan_error());
|
|
}
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
}
|
|
|
|
#ifdef FE_INVALID
|
|
if (feclearexcept(FE_INVALID)) throw fenv_error();
|
|
#endif
|
|
|
|
for (t_index j=0; j<N-1; ++j) {
|
|
if (NN_chain_tip <= 3) {
|
|
NN_chain[0] = idx1 = active_nodes.start;
|
|
NN_chain_tip = 1;
|
|
|
|
idx2 = active_nodes.succ[idx1];
|
|
min = D_(idx1,idx2);
|
|
for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i]) {
|
|
if (D_(idx1,i) < min) {
|
|
min = D_(idx1,i);
|
|
idx2 = i;
|
|
}
|
|
}
|
|
} // a: idx1 b: idx2
|
|
else {
|
|
NN_chain_tip -= 3;
|
|
idx1 = NN_chain[NN_chain_tip-1];
|
|
idx2 = NN_chain[NN_chain_tip];
|
|
min = idx1<idx2 ? D_(idx1,idx2) : D_(idx2,idx1);
|
|
} // a: idx1 b: idx2
|
|
|
|
do {
|
|
NN_chain[NN_chain_tip] = idx2;
|
|
|
|
for (i=active_nodes.start; i<idx2; i=active_nodes.succ[i]) {
|
|
if (D_(i,idx2) < min) {
|
|
min = D_(i,idx2);
|
|
idx1 = i;
|
|
}
|
|
}
|
|
for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i]) {
|
|
if (D_(idx2,i) < min) {
|
|
min = D_(idx2,i);
|
|
idx1 = i;
|
|
}
|
|
}
|
|
|
|
idx2 = idx1;
|
|
idx1 = NN_chain[NN_chain_tip++];
|
|
|
|
} while (idx2 != NN_chain[NN_chain_tip-2]);
|
|
|
|
Z2.append(idx1, idx2, min);
|
|
|
|
if (idx1>idx2) {
|
|
t_index tmp = idx1;
|
|
idx1 = idx2;
|
|
idx2 = tmp;
|
|
}
|
|
|
|
if (method==METHOD_METR_AVERAGE ||
|
|
method==METHOD_METR_WARD) {
|
|
size1 = static_cast<t_float>(members[idx1]);
|
|
size2 = static_cast<t_float>(members[idx2]);
|
|
members[idx2] += members[idx1];
|
|
}
|
|
|
|
// Remove the smaller index from the valid indices (active_nodes).
|
|
active_nodes.remove(idx1);
|
|
|
|
switch (method) {
|
|
case METHOD_METR_SINGLE:
|
|
/*
|
|
Single linkage.
|
|
|
|
Characteristic: new distances are never longer than the old distances.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (i=active_nodes.start; i<idx1; i=active_nodes.succ[i])
|
|
f_single(&D_(i, idx2), D_(i, idx1) );
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; i<idx2; i=active_nodes.succ[i])
|
|
f_single(&D_(i, idx2), D_(idx1, i) );
|
|
// Update the distance matrix in the range (idx2, N).
|
|
for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i])
|
|
f_single(&D_(idx2, i), D_(idx1, i) );
|
|
break;
|
|
|
|
case METHOD_METR_COMPLETE:
|
|
/*
|
|
Complete linkage.
|
|
|
|
Characteristic: new distances are never shorter than the old distances.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (i=active_nodes.start; i<idx1; i=active_nodes.succ[i])
|
|
f_complete(&D_(i, idx2), D_(i, idx1) );
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; i<idx2; i=active_nodes.succ[i])
|
|
f_complete(&D_(i, idx2), D_(idx1, i) );
|
|
// Update the distance matrix in the range (idx2, N).
|
|
for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i])
|
|
f_complete(&D_(idx2, i), D_(idx1, i) );
|
|
break;
|
|
|
|
case METHOD_METR_AVERAGE: {
|
|
/*
|
|
Average linkage.
|
|
|
|
Shorter and longer distances can occur.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
t_float s = size1/(size1+size2);
|
|
t_float t = size2/(size1+size2);
|
|
for (i=active_nodes.start; i<idx1; i=active_nodes.succ[i])
|
|
f_average(&D_(i, idx2), D_(i, idx1), s, t );
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; i<idx2; i=active_nodes.succ[i])
|
|
f_average(&D_(i, idx2), D_(idx1, i), s, t );
|
|
// Update the distance matrix in the range (idx2, N).
|
|
for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i])
|
|
f_average(&D_(idx2, i), D_(idx1, i), s, t );
|
|
break;
|
|
}
|
|
|
|
case METHOD_METR_WEIGHTED:
|
|
/*
|
|
Weighted linkage.
|
|
|
|
Shorter and longer distances can occur.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (i=active_nodes.start; i<idx1; i=active_nodes.succ[i])
|
|
f_weighted(&D_(i, idx2), D_(i, idx1) );
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; i<idx2; i=active_nodes.succ[i])
|
|
f_weighted(&D_(i, idx2), D_(idx1, i) );
|
|
// Update the distance matrix in the range (idx2, N).
|
|
for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i])
|
|
f_weighted(&D_(idx2, i), D_(idx1, i) );
|
|
break;
|
|
|
|
case METHOD_METR_WARD:
|
|
/*
|
|
Ward linkage.
|
|
|
|
Shorter and longer distances can occur, not smaller than min(d1,d2)
|
|
but maybe bigger than max(d1,d2).
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
//t_float v = static_cast<t_float>(members[i]);
|
|
for (i=active_nodes.start; i<idx1; i=active_nodes.succ[i])
|
|
f_ward(&D_(i, idx2), D_(i, idx1), min,
|
|
size1, size2, static_cast<t_float>(members[i]) );
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; i<idx2; i=active_nodes.succ[i])
|
|
f_ward(&D_(i, idx2), D_(idx1, i), min,
|
|
size1, size2, static_cast<t_float>(members[i]) );
|
|
// Update the distance matrix in the range (idx2, N).
|
|
for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i])
|
|
f_ward(&D_(idx2, i), D_(idx1, i), min,
|
|
size1, size2, static_cast<t_float>(members[i]) );
|
|
break;
|
|
|
|
default:
|
|
throw std::runtime_error(std::string("Invalid method."));
|
|
}
|
|
}
|
|
#ifdef FE_INVALID
|
|
if (fetestexcept(FE_INVALID)) throw fenv_error();
|
|
#endif
|
|
}
|
|
|
|
class binary_min_heap {
|
|
/*
|
|
Class for a binary min-heap. The data resides in an array A. The elements of
|
|
A are not changed but two lists I and R of indices are generated which point
|
|
to elements of A and backwards.
|
|
|
|
The heap tree structure is
|
|
|
|
H[2*i+1] H[2*i+2]
|
|
\ /
|
|
\ /
|
|
≤ ≤
|
|
\ /
|
|
\ /
|
|
H[i]
|
|
|
|
where the children must be less or equal than their parent. Thus, H[0]
|
|
contains the minimum. The lists I and R are made such that H[i] = A[I[i]]
|
|
and R[I[i]] = i.
|
|
|
|
This implementation is not designed to handle NaN values.
|
|
*/
|
|
private:
|
|
t_float * const A;
|
|
t_index size;
|
|
auto_array_ptr<t_index> I;
|
|
auto_array_ptr<t_index> R;
|
|
|
|
// no default constructor
|
|
binary_min_heap();
|
|
// noncopyable
|
|
binary_min_heap(binary_min_heap const &);
|
|
binary_min_heap & operator=(binary_min_heap const &);
|
|
|
|
public:
|
|
binary_min_heap(t_float * const A_, const t_index size_)
|
|
: A(A_), size(size_), I(size), R(size)
|
|
{ // Allocate memory and initialize the lists I and R to the identity. This
|
|
// does not make it a heap. Call heapify afterwards!
|
|
for (t_index i=0; i<size; ++i)
|
|
R[i] = I[i] = i;
|
|
}
|
|
|
|
binary_min_heap(t_float * const A_, const t_index size1, const t_index size2,
|
|
const t_index start)
|
|
: A(A_), size(size1), I(size1), R(size2)
|
|
{ // Allocate memory and initialize the lists I and R to the identity. This
|
|
// does not make it a heap. Call heapify afterwards!
|
|
for (t_index i=0; i<size; ++i) {
|
|
R[i+start] = i;
|
|
I[i] = i + start;
|
|
}
|
|
}
|
|
|
|
~binary_min_heap() {}
|
|
|
|
void heapify() {
|
|
// Arrange the indices I and R so that H[i] := A[I[i]] satisfies the heap
|
|
// condition H[i] < H[2*i+1] and H[i] < H[2*i+2] for each i.
|
|
//
|
|
// Complexity: Θ(size)
|
|
// Reference: Cormen, Leiserson, Rivest, Stein, Introduction to Algorithms,
|
|
// 3rd ed., 2009, Section 6.3 “Building a heap”
|
|
t_index idx;
|
|
for (idx=(size>>1); idx>0; ) {
|
|
--idx;
|
|
update_geq_(idx);
|
|
}
|
|
}
|
|
|
|
inline t_index argmin() const {
|
|
// Return the minimal element.
|
|
return I[0];
|
|
}
|
|
|
|
void heap_pop() {
|
|
// Remove the minimal element from the heap.
|
|
--size;
|
|
I[0] = I[size];
|
|
R[I[0]] = 0;
|
|
update_geq_(0);
|
|
}
|
|
|
|
void remove(t_index idx) {
|
|
// Remove an element from the heap.
|
|
--size;
|
|
R[I[size]] = R[idx];
|
|
I[R[idx]] = I[size];
|
|
if ( H(size)<=A[idx] ) {
|
|
update_leq_(R[idx]);
|
|
}
|
|
else {
|
|
update_geq_(R[idx]);
|
|
}
|
|
}
|
|
|
|
void replace ( const t_index idxold, const t_index idxnew,
|
|
const t_float val) {
|
|
R[idxnew] = R[idxold];
|
|
I[R[idxnew]] = idxnew;
|
|
if (val<=A[idxold])
|
|
update_leq(idxnew, val);
|
|
else
|
|
update_geq(idxnew, val);
|
|
}
|
|
|
|
void update ( const t_index idx, const t_float val ) const {
|
|
// Update the element A[i] with val and re-arrange the indices to preserve
|
|
// the heap condition.
|
|
if (val<=A[idx])
|
|
update_leq(idx, val);
|
|
else
|
|
update_geq(idx, val);
|
|
}
|
|
|
|
void update_leq ( const t_index idx, const t_float val ) const {
|
|
// Use this when the new value is not more than the old value.
|
|
A[idx] = val;
|
|
update_leq_(R[idx]);
|
|
}
|
|
|
|
void update_geq ( const t_index idx, const t_float val ) const {
|
|
// Use this when the new value is not less than the old value.
|
|
A[idx] = val;
|
|
update_geq_(R[idx]);
|
|
}
|
|
|
|
private:
|
|
void update_leq_ (t_index i) const {
|
|
t_index j;
|
|
for ( ; (i>0) && ( H(i)<H(j=(i-1)>>1) ); i=j)
|
|
heap_swap(i,j);
|
|
}
|
|
|
|
void update_geq_ (t_index i) const {
|
|
t_index j;
|
|
for ( ; (j=2*i+1)<size; i=j) {
|
|
if ( H(j)>=H(i) ) {
|
|
++j;
|
|
if ( j>=size || H(j)>=H(i) ) break;
|
|
}
|
|
else if ( j+1<size && H(j+1)<H(j) ) ++j;
|
|
heap_swap(i, j);
|
|
}
|
|
}
|
|
|
|
void heap_swap(const t_index i, const t_index j) const {
|
|
// Swap two indices.
|
|
t_index tmp = I[i];
|
|
I[i] = I[j];
|
|
I[j] = tmp;
|
|
R[I[i]] = i;
|
|
R[I[j]] = j;
|
|
}
|
|
|
|
inline t_float H(const t_index i) const {
|
|
return A[I[i]];
|
|
}
|
|
|
|
};
|
|
|
|
template <method_codes method, typename t_members>
|
|
static void generic_linkage(const t_index N, t_float * const D, t_members * const members, cluster_result & Z2) {
|
|
/*
|
|
N: integer, number of data points
|
|
D: condensed distance matrix N*(N-1)/2
|
|
Z2: output data structure
|
|
*/
|
|
|
|
const t_index N_1 = N-1;
|
|
t_index i, j; // loop variables
|
|
t_index idx1, idx2; // row and column indices
|
|
|
|
auto_array_ptr<t_index> n_nghbr(N_1); // array of nearest neighbors
|
|
auto_array_ptr<t_float> mindist(N_1); // distances to the nearest neighbors
|
|
auto_array_ptr<t_index> row_repr(N); // row_repr[i]: node number that the
|
|
// i-th row represents
|
|
doubly_linked_list active_nodes(N);
|
|
binary_min_heap nn_distances(&*mindist, N_1); // minimum heap structure for
|
|
// the distance to the nearest neighbor of each point
|
|
t_index node1, node2; // node numbers in the output
|
|
t_float size1, size2; // and their cardinalities
|
|
|
|
t_float min; // minimum and row index for nearest-neighbor search
|
|
t_index idx;
|
|
|
|
for (i=0; i<N; ++i)
|
|
// Build a list of row ↔ node label assignments.
|
|
// Initially i ↦ i
|
|
row_repr[i] = i;
|
|
|
|
// Initialize the minimal distances:
|
|
// Find the nearest neighbor of each point.
|
|
// n_nghbr[i] = argmin_{j>i} D(i,j) for i in range(N-1)
|
|
t_float const * DD = D;
|
|
for (i=0; i<N_1; ++i) {
|
|
min = std::numeric_limits<t_float>::infinity();
|
|
for (idx=j=i+1; j<N; ++j, ++DD) {
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (*DD<min) {
|
|
min = *DD;
|
|
idx = j;
|
|
}
|
|
else if (fc_isnan(*DD))
|
|
throw(nan_error());
|
|
}
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
mindist[i] = min;
|
|
n_nghbr[i] = idx;
|
|
}
|
|
|
|
// Put the minimal distances into a heap structure to make the repeated
|
|
// global minimum searches fast.
|
|
nn_distances.heapify();
|
|
|
|
#ifdef FE_INVALID
|
|
if (feclearexcept(FE_INVALID)) throw fenv_error();
|
|
#endif
|
|
|
|
// Main loop: We have N-1 merging steps.
|
|
for (i=0; i<N_1; ++i) {
|
|
/*
|
|
Here is a special feature that allows fast bookkeeping and updates of the
|
|
minimal distances.
|
|
|
|
mindist[i] stores a lower bound on the minimum distance of the point i to
|
|
all points of higher index:
|
|
|
|
mindist[i] ≥ min_{j>i} D(i,j)
|
|
|
|
Normally, we have equality. However, this minimum may become invalid due
|
|
to the updates in the distance matrix. The rules are:
|
|
|
|
1) If mindist[i] is equal to D(i, n_nghbr[i]), this is the correct
|
|
minimum and n_nghbr[i] is a nearest neighbor.
|
|
|
|
2) If mindist[i] is smaller than D(i, n_nghbr[i]), this might not be the
|
|
correct minimum. The minimum needs to be recomputed.
|
|
|
|
3) mindist[i] is never bigger than the true minimum. Hence, we never
|
|
miss the true minimum if we take the smallest mindist entry,
|
|
re-compute the value if necessary (thus maybe increasing it) and
|
|
looking for the now smallest mindist entry until a valid minimal
|
|
entry is found. This step is done in the lines below.
|
|
|
|
The update process for D below takes care that these rules are
|
|
fulfilled. This makes sure that the minima in the rows D(i,i+1:)of D are
|
|
re-calculated when necessary but re-calculation is avoided whenever
|
|
possible.
|
|
|
|
The re-calculation of the minima makes the worst-case runtime of this
|
|
algorithm cubic in N. We avoid this whenever possible, and in most cases
|
|
the runtime appears to be quadratic.
|
|
*/
|
|
idx1 = nn_distances.argmin();
|
|
if (method != METHOD_METR_SINGLE) {
|
|
while ( mindist[idx1] < D_(idx1, n_nghbr[idx1]) ) {
|
|
// Recompute the minimum mindist[idx1] and n_nghbr[idx1].
|
|
n_nghbr[idx1] = j = active_nodes.succ[idx1]; // exists, maximally N-1
|
|
min = D_(idx1,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
if (D_(idx1,j)<min) {
|
|
min = D_(idx1,j);
|
|
n_nghbr[idx1] = j;
|
|
}
|
|
}
|
|
/* Update the heap with the new true minimum and search for the
|
|
(possibly different) minimal entry. */
|
|
nn_distances.update_geq(idx1, min);
|
|
idx1 = nn_distances.argmin();
|
|
}
|
|
}
|
|
|
|
nn_distances.heap_pop(); // Remove the current minimum from the heap.
|
|
idx2 = n_nghbr[idx1];
|
|
|
|
// Write the newly found minimal pair of nodes to the output array.
|
|
node1 = row_repr[idx1];
|
|
node2 = row_repr[idx2];
|
|
|
|
if (method==METHOD_METR_AVERAGE ||
|
|
method==METHOD_METR_WARD ||
|
|
method==METHOD_METR_CENTROID) {
|
|
size1 = static_cast<t_float>(members[idx1]);
|
|
size2 = static_cast<t_float>(members[idx2]);
|
|
members[idx2] += members[idx1];
|
|
}
|
|
Z2.append(node1, node2, mindist[idx1]);
|
|
|
|
// Remove idx1 from the list of active indices (active_nodes).
|
|
active_nodes.remove(idx1);
|
|
// Index idx2 now represents the new (merged) node with label N+i.
|
|
row_repr[idx2] = N+i;
|
|
|
|
// Update the distance matrix
|
|
switch (method) {
|
|
case METHOD_METR_SINGLE:
|
|
/*
|
|
Single linkage.
|
|
|
|
Characteristic: new distances are never longer than the old distances.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
f_single(&D_(j, idx2), D_(j, idx1));
|
|
if (n_nghbr[j] == idx1)
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; j<idx2; j=active_nodes.succ[j]) {
|
|
f_single(&D_(j, idx2), D_(idx1, j));
|
|
// If the new value is below the old minimum in a row, update
|
|
// the mindist and n_nghbr arrays.
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
}
|
|
// Update the distance matrix in the range (idx2, N).
|
|
// Recompute the minimum mindist[idx2] and n_nghbr[idx2].
|
|
if (idx2<N_1) {
|
|
min = mindist[idx2];
|
|
for (j=active_nodes.succ[idx2]; j<N; j=active_nodes.succ[j]) {
|
|
f_single(&D_(idx2, j), D_(idx1, j) );
|
|
if (D_(idx2, j) < min) {
|
|
n_nghbr[idx2] = j;
|
|
min = D_(idx2, j);
|
|
}
|
|
}
|
|
nn_distances.update_leq(idx2, min);
|
|
}
|
|
break;
|
|
|
|
case METHOD_METR_COMPLETE:
|
|
/*
|
|
Complete linkage.
|
|
|
|
Characteristic: new distances are never shorter than the old distances.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
f_complete(&D_(j, idx2), D_(j, idx1) );
|
|
if (n_nghbr[j] == idx1)
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; j<idx2; j=active_nodes.succ[j])
|
|
f_complete(&D_(j, idx2), D_(idx1, j) );
|
|
// Update the distance matrix in the range (idx2, N).
|
|
for (j=active_nodes.succ[idx2]; j<N; j=active_nodes.succ[j])
|
|
f_complete(&D_(idx2, j), D_(idx1, j) );
|
|
break;
|
|
|
|
case METHOD_METR_AVERAGE: {
|
|
/*
|
|
Average linkage.
|
|
|
|
Shorter and longer distances can occur.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
t_float s = size1/(size1+size2);
|
|
t_float t = size2/(size1+size2);
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
f_average(&D_(j, idx2), D_(j, idx1), s, t);
|
|
if (n_nghbr[j] == idx1)
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; j<idx2; j=active_nodes.succ[j]) {
|
|
f_average(&D_(j, idx2), D_(idx1, j), s, t);
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
}
|
|
// Update the distance matrix in the range (idx2, N).
|
|
if (idx2<N_1) {
|
|
n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
|
|
f_average(&D_(idx2, j), D_(idx1, j), s, t);
|
|
min = D_(idx2,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
f_average(&D_(idx2, j), D_(idx1, j), s, t);
|
|
if (D_(idx2,j) < min) {
|
|
min = D_(idx2,j);
|
|
n_nghbr[idx2] = j;
|
|
}
|
|
}
|
|
nn_distances.update(idx2, min);
|
|
}
|
|
break;
|
|
}
|
|
|
|
case METHOD_METR_WEIGHTED:
|
|
/*
|
|
Weighted linkage.
|
|
|
|
Shorter and longer distances can occur.
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
f_weighted(&D_(j, idx2), D_(j, idx1) );
|
|
if (n_nghbr[j] == idx1)
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; j<idx2; j=active_nodes.succ[j]) {
|
|
f_weighted(&D_(j, idx2), D_(idx1, j) );
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
}
|
|
// Update the distance matrix in the range (idx2, N).
|
|
if (idx2<N_1) {
|
|
n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
|
|
f_weighted(&D_(idx2, j), D_(idx1, j) );
|
|
min = D_(idx2,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
f_weighted(&D_(idx2, j), D_(idx1, j) );
|
|
if (D_(idx2,j) < min) {
|
|
min = D_(idx2,j);
|
|
n_nghbr[idx2] = j;
|
|
}
|
|
}
|
|
nn_distances.update(idx2, min);
|
|
}
|
|
break;
|
|
|
|
case METHOD_METR_WARD:
|
|
/*
|
|
Ward linkage.
|
|
|
|
Shorter and longer distances can occur, not smaller than min(d1,d2)
|
|
but maybe bigger than max(d1,d2).
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
f_ward(&D_(j, idx2), D_(j, idx1), mindist[idx1],
|
|
size1, size2, static_cast<t_float>(members[j]) );
|
|
if (n_nghbr[j] == idx1)
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; j<idx2; j=active_nodes.succ[j]) {
|
|
f_ward(&D_(j, idx2), D_(idx1, j), mindist[idx1], size1, size2,
|
|
static_cast<t_float>(members[j]) );
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
}
|
|
// Update the distance matrix in the range (idx2, N).
|
|
if (idx2<N_1) {
|
|
n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
|
|
f_ward(&D_(idx2, j), D_(idx1, j), mindist[idx1],
|
|
size1, size2, static_cast<t_float>(members[j]) );
|
|
min = D_(idx2,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
f_ward(&D_(idx2, j), D_(idx1, j), mindist[idx1],
|
|
size1, size2, static_cast<t_float>(members[j]) );
|
|
if (D_(idx2,j) < min) {
|
|
min = D_(idx2,j);
|
|
n_nghbr[idx2] = j;
|
|
}
|
|
}
|
|
nn_distances.update(idx2, min);
|
|
}
|
|
break;
|
|
|
|
case METHOD_METR_CENTROID: {
|
|
/*
|
|
Centroid linkage.
|
|
|
|
Shorter and longer distances can occur, not bigger than max(d1,d2)
|
|
but maybe smaller than min(d1,d2).
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
t_float s = size1/(size1+size2);
|
|
t_float t = size2/(size1+size2);
|
|
t_float stc = s*t*mindist[idx1];
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
f_centroid(&D_(j, idx2), D_(j, idx1), stc, s, t);
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
else if (n_nghbr[j] == idx1)
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; j<idx2; j=active_nodes.succ[j]) {
|
|
f_centroid(&D_(j, idx2), D_(idx1, j), stc, s, t);
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
}
|
|
// Update the distance matrix in the range (idx2, N).
|
|
if (idx2<N_1) {
|
|
n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
|
|
f_centroid(&D_(idx2, j), D_(idx1, j), stc, s, t);
|
|
min = D_(idx2,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
f_centroid(&D_(idx2, j), D_(idx1, j), stc, s, t);
|
|
if (D_(idx2,j) < min) {
|
|
min = D_(idx2,j);
|
|
n_nghbr[idx2] = j;
|
|
}
|
|
}
|
|
nn_distances.update(idx2, min);
|
|
}
|
|
break;
|
|
}
|
|
|
|
case METHOD_METR_MEDIAN: {
|
|
/*
|
|
Median linkage.
|
|
|
|
Shorter and longer distances can occur, not bigger than max(d1,d2)
|
|
but maybe smaller than min(d1,d2).
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
t_float c_4 = mindist[idx1]*.25;
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
f_median(&D_(j, idx2), D_(j, idx1), c_4 );
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
else if (n_nghbr[j] == idx1)
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for (; j<idx2; j=active_nodes.succ[j]) {
|
|
f_median(&D_(j, idx2), D_(idx1, j), c_4 );
|
|
if (D_(j, idx2) < mindist[j]) {
|
|
nn_distances.update_leq(j, D_(j, idx2));
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
}
|
|
// Update the distance matrix in the range (idx2, N).
|
|
if (idx2<N_1) {
|
|
n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
|
|
f_median(&D_(idx2, j), D_(idx1, j), c_4 );
|
|
min = D_(idx2,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
f_median(&D_(idx2, j), D_(idx1, j), c_4 );
|
|
if (D_(idx2,j) < min) {
|
|
min = D_(idx2,j);
|
|
n_nghbr[idx2] = j;
|
|
}
|
|
}
|
|
nn_distances.update(idx2, min);
|
|
}
|
|
break;
|
|
}
|
|
|
|
default:
|
|
throw std::runtime_error(std::string("Invalid method."));
|
|
}
|
|
}
|
|
#ifdef FE_INVALID
|
|
if (fetestexcept(FE_INVALID)) throw fenv_error();
|
|
#endif
|
|
}
|
|
|
|
/*
|
|
Clustering methods for vector data
|
|
*/
|
|
|
|
template <typename t_dissimilarity>
|
|
static void MST_linkage_core_vector(const t_index N,
|
|
t_dissimilarity & dist,
|
|
cluster_result & Z2) {
|
|
/*
|
|
N: integer, number of data points
|
|
dist: function pointer to the metric
|
|
Z2: output data structure
|
|
|
|
The basis of this algorithm is an algorithm by Rohlf:
|
|
|
|
F. James Rohlf, Hierarchical clustering using the minimum spanning tree,
|
|
The Computer Journal, vol. 16, 1973, p. 93–95.
|
|
*/
|
|
t_index i;
|
|
t_index idx2;
|
|
doubly_linked_list active_nodes(N);
|
|
auto_array_ptr<t_float> d(N);
|
|
|
|
t_index prev_node;
|
|
t_float min;
|
|
|
|
// first iteration
|
|
idx2 = 1;
|
|
min = std::numeric_limits<t_float>::infinity();
|
|
for (i=1; i<N; ++i) {
|
|
d[i] = dist(0,i);
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (d[i] < min) {
|
|
min = d[i];
|
|
idx2 = i;
|
|
}
|
|
else if (fc_isnan(d[i]))
|
|
throw (nan_error());
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
}
|
|
|
|
Z2.append(0, idx2, min);
|
|
|
|
for (t_index j=1; j<N-1; ++j) {
|
|
prev_node = idx2;
|
|
active_nodes.remove(prev_node);
|
|
|
|
idx2 = active_nodes.succ[0];
|
|
min = d[idx2];
|
|
|
|
for (i=idx2; i<N; i=active_nodes.succ[i]) {
|
|
t_float tmp = dist(i, prev_node);
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
if (d[i] > tmp)
|
|
d[i] = tmp;
|
|
else if (fc_isnan(tmp))
|
|
throw (nan_error());
|
|
#if HAVE_DIAGNOSTIC
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
if (d[i] < min) {
|
|
min = d[i];
|
|
idx2 = i;
|
|
}
|
|
}
|
|
Z2.append(prev_node, idx2, min);
|
|
}
|
|
}
|
|
|
|
template <method_codes_vector method, typename t_dissimilarity>
|
|
static void generic_linkage_vector(const t_index N,
|
|
t_dissimilarity & dist,
|
|
cluster_result & Z2) {
|
|
/*
|
|
N: integer, number of data points
|
|
dist: function pointer to the metric
|
|
Z2: output data structure
|
|
|
|
This algorithm is valid for the distance update methods
|
|
"Ward", "centroid" and "median" only!
|
|
*/
|
|
const t_index N_1 = N-1;
|
|
t_index i, j; // loop variables
|
|
t_index idx1, idx2; // row and column indices
|
|
|
|
auto_array_ptr<t_index> n_nghbr(N_1); // array of nearest neighbors
|
|
auto_array_ptr<t_float> mindist(N_1); // distances to the nearest neighbors
|
|
auto_array_ptr<t_index> row_repr(N); // row_repr[i]: node number that the
|
|
// i-th row represents
|
|
doubly_linked_list active_nodes(N);
|
|
binary_min_heap nn_distances(&*mindist, N_1); // minimum heap structure for
|
|
// the distance to the nearest neighbor of each point
|
|
t_index node1, node2; // node numbers in the output
|
|
t_float min; // minimum and row index for nearest-neighbor search
|
|
|
|
for (i=0; i<N; ++i)
|
|
// Build a list of row ↔ node label assignments.
|
|
// Initially i ↦ i
|
|
row_repr[i] = i;
|
|
|
|
// Initialize the minimal distances:
|
|
// Find the nearest neighbor of each point.
|
|
// n_nghbr[i] = argmin_{j>i} D(i,j) for i in range(N-1)
|
|
for (i=0; i<N_1; ++i) {
|
|
min = std::numeric_limits<t_float>::infinity();
|
|
t_index idx;
|
|
for (idx=j=i+1; j<N; ++j) {
|
|
t_float tmp;
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
tmp = dist.ward_initial(i,j);
|
|
break;
|
|
default:
|
|
tmp = dist.template sqeuclidean<true>(i,j);
|
|
}
|
|
if (tmp<min) {
|
|
min = tmp;
|
|
idx = j;
|
|
}
|
|
}
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
mindist[i] = t_dissimilarity::ward_initial_conversion(min);
|
|
break;
|
|
default:
|
|
mindist[i] = min;
|
|
}
|
|
n_nghbr[i] = idx;
|
|
}
|
|
|
|
// Put the minimal distances into a heap structure to make the repeated
|
|
// global minimum searches fast.
|
|
nn_distances.heapify();
|
|
|
|
// Main loop: We have N-1 merging steps.
|
|
for (i=0; i<N_1; ++i) {
|
|
idx1 = nn_distances.argmin();
|
|
|
|
while ( active_nodes.is_inactive(n_nghbr[idx1]) ) {
|
|
// Recompute the minimum mindist[idx1] and n_nghbr[idx1].
|
|
n_nghbr[idx1] = j = active_nodes.succ[idx1]; // exists, maximally N-1
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
min = dist.ward(idx1,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
t_float const tmp = dist.ward(idx1,j);
|
|
if (tmp<min) {
|
|
min = tmp;
|
|
n_nghbr[idx1] = j;
|
|
}
|
|
}
|
|
break;
|
|
default:
|
|
min = dist.template sqeuclidean<true>(idx1,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
t_float const tmp = dist.template sqeuclidean<true>(idx1,j);
|
|
if (tmp<min) {
|
|
min = tmp;
|
|
n_nghbr[idx1] = j;
|
|
}
|
|
}
|
|
}
|
|
/* Update the heap with the new true minimum and search for the (possibly
|
|
different) minimal entry. */
|
|
nn_distances.update_geq(idx1, min);
|
|
idx1 = nn_distances.argmin();
|
|
}
|
|
|
|
nn_distances.heap_pop(); // Remove the current minimum from the heap.
|
|
idx2 = n_nghbr[idx1];
|
|
|
|
// Write the newly found minimal pair of nodes to the output array.
|
|
node1 = row_repr[idx1];
|
|
node2 = row_repr[idx2];
|
|
|
|
Z2.append(node1, node2, mindist[idx1]);
|
|
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
case METHOD_VECTOR_CENTROID:
|
|
dist.merge_inplace(idx1, idx2);
|
|
break;
|
|
case METHOD_VECTOR_MEDIAN:
|
|
dist.merge_inplace_weighted(idx1, idx2);
|
|
break;
|
|
default:
|
|
throw std::runtime_error(std::string("Invalid method."));
|
|
}
|
|
|
|
// Index idx2 now represents the new (merged) node with label N+i.
|
|
row_repr[idx2] = N+i;
|
|
// Remove idx1 from the list of active indices (active_nodes).
|
|
active_nodes.remove(idx1); // TBD later!!!
|
|
|
|
// Update the distance matrix
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
/*
|
|
Ward linkage.
|
|
|
|
Shorter and longer distances can occur, not smaller than min(d1,d2)
|
|
but maybe bigger than max(d1,d2).
|
|
*/
|
|
// Update the distance matrix in the range [start, idx1).
|
|
for (j=active_nodes.start; j<idx1; j=active_nodes.succ[j]) {
|
|
if (n_nghbr[j] == idx2) {
|
|
n_nghbr[j] = idx1; // invalidate
|
|
}
|
|
}
|
|
// Update the distance matrix in the range (idx1, idx2).
|
|
for ( ; j<idx2; j=active_nodes.succ[j]) {
|
|
t_float const tmp = dist.ward(j, idx2);
|
|
if (tmp < mindist[j]) {
|
|
nn_distances.update_leq(j, tmp);
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
else if (n_nghbr[j]==idx2) {
|
|
n_nghbr[j] = idx1; // invalidate
|
|
}
|
|
}
|
|
// Find the nearest neighbor for idx2.
|
|
if (idx2<N_1) {
|
|
n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
|
|
min = dist.ward(idx2,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
t_float const tmp = dist.ward(idx2,j);
|
|
if (tmp < min) {
|
|
min = tmp;
|
|
n_nghbr[idx2] = j;
|
|
}
|
|
}
|
|
nn_distances.update(idx2, min);
|
|
}
|
|
break;
|
|
|
|
default:
|
|
/*
|
|
Centroid and median linkage.
|
|
|
|
Shorter and longer distances can occur, not bigger than max(d1,d2)
|
|
but maybe smaller than min(d1,d2).
|
|
*/
|
|
for (j=active_nodes.start; j<idx2; j=active_nodes.succ[j]) {
|
|
t_float const tmp = dist.template sqeuclidean<true>(j, idx2);
|
|
if (tmp < mindist[j]) {
|
|
nn_distances.update_leq(j, tmp);
|
|
n_nghbr[j] = idx2;
|
|
}
|
|
else if (n_nghbr[j] == idx2)
|
|
n_nghbr[j] = idx1; // invalidate
|
|
}
|
|
// Find the nearest neighbor for idx2.
|
|
if (idx2<N_1) {
|
|
n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
|
|
min = dist.template sqeuclidean<true>(idx2,j);
|
|
for (j=active_nodes.succ[j]; j<N; j=active_nodes.succ[j]) {
|
|
t_float const tmp = dist.template sqeuclidean<true>(idx2, j);
|
|
if (tmp < min) {
|
|
min = tmp;
|
|
n_nghbr[idx2] = j;
|
|
}
|
|
}
|
|
nn_distances.update(idx2, min);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template <method_codes_vector method, typename t_dissimilarity>
|
|
static void generic_linkage_vector_alternative(const t_index N,
|
|
t_dissimilarity & dist,
|
|
cluster_result & Z2) {
|
|
/*
|
|
N: integer, number of data points
|
|
dist: function pointer to the metric
|
|
Z2: output data structure
|
|
|
|
This algorithm is valid for the distance update methods
|
|
"Ward", "centroid" and "median" only!
|
|
*/
|
|
const t_index N_1 = N-1;
|
|
t_index i, j=0; // loop variables
|
|
t_index idx1, idx2; // row and column indices
|
|
|
|
auto_array_ptr<t_index> n_nghbr(2*N-2); // array of nearest neighbors
|
|
auto_array_ptr<t_float> mindist(2*N-2); // distances to the nearest neighbors
|
|
|
|
doubly_linked_list active_nodes(N+N_1);
|
|
binary_min_heap nn_distances(&*mindist, N_1, 2*N-2, 1); // minimum heap
|
|
// structure for the distance to the nearest neighbor of each point
|
|
|
|
t_float min; // minimum for nearest-neighbor searches
|
|
|
|
// Initialize the minimal distances:
|
|
// Find the nearest neighbor of each point.
|
|
// n_nghbr[i] = argmin_{j>i} D(i,j) for i in range(N-1)
|
|
for (i=1; i<N; ++i) {
|
|
min = std::numeric_limits<t_float>::infinity();
|
|
t_index idx;
|
|
for (idx=j=0; j<i; ++j) {
|
|
t_float tmp;
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
tmp = dist.ward_initial(i,j);
|
|
break;
|
|
default:
|
|
tmp = dist.template sqeuclidean<true>(i,j);
|
|
}
|
|
if (tmp<min) {
|
|
min = tmp;
|
|
idx = j;
|
|
}
|
|
}
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
mindist[i] = t_dissimilarity::ward_initial_conversion(min);
|
|
break;
|
|
default:
|
|
mindist[i] = min;
|
|
}
|
|
n_nghbr[i] = idx;
|
|
}
|
|
|
|
// Put the minimal distances into a heap structure to make the repeated
|
|
// global minimum searches fast.
|
|
nn_distances.heapify();
|
|
|
|
// Main loop: We have N-1 merging steps.
|
|
for (i=N; i<N+N_1; ++i) {
|
|
/*
|
|
The bookkeeping is different from the "stored matrix approach" algorithm
|
|
generic_linkage.
|
|
|
|
mindist[i] stores a lower bound on the minimum distance of the point i to
|
|
all points of *lower* index:
|
|
|
|
mindist[i] ≥ min_{j<i} D(i,j)
|
|
|
|
Moreover, new nodes do not re-use one of the old indices, but they are
|
|
given a new, unique index (SciPy convention: initial nodes are 0,…,N−1,
|
|
new nodes are N,…,2N−2).
|
|
|
|
Invalid nearest neighbors are not recognized by the fact that the stored
|
|
distance is smaller than the actual distance, but the list active_nodes
|
|
maintains a flag whether a node is inactive. If n_nghbr[i] points to an
|
|
active node, the entries nn_distances[i] and n_nghbr[i] are valid,
|
|
otherwise they must be recomputed.
|
|
*/
|
|
idx1 = nn_distances.argmin();
|
|
while ( active_nodes.is_inactive(n_nghbr[idx1]) ) {
|
|
// Recompute the minimum mindist[idx1] and n_nghbr[idx1].
|
|
n_nghbr[idx1] = j = active_nodes.start;
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
min = dist.ward_extended(idx1,j);
|
|
for (j=active_nodes.succ[j]; j<idx1; j=active_nodes.succ[j]) {
|
|
t_float tmp = dist.ward_extended(idx1,j);
|
|
if (tmp<min) {
|
|
min = tmp;
|
|
n_nghbr[idx1] = j;
|
|
}
|
|
}
|
|
break;
|
|
default:
|
|
min = dist.sqeuclidean_extended(idx1,j);
|
|
for (j=active_nodes.succ[j]; j<idx1; j=active_nodes.succ[j]) {
|
|
t_float const tmp = dist.sqeuclidean_extended(idx1,j);
|
|
if (tmp<min) {
|
|
min = tmp;
|
|
n_nghbr[idx1] = j;
|
|
}
|
|
}
|
|
}
|
|
/* Update the heap with the new true minimum and search for the (possibly
|
|
different) minimal entry. */
|
|
nn_distances.update_geq(idx1, min);
|
|
idx1 = nn_distances.argmin();
|
|
}
|
|
|
|
idx2 = n_nghbr[idx1];
|
|
active_nodes.remove(idx1);
|
|
active_nodes.remove(idx2);
|
|
|
|
Z2.append(idx1, idx2, mindist[idx1]);
|
|
|
|
if (i<2*N_1) {
|
|
switch (method) {
|
|
case METHOD_VECTOR_WARD:
|
|
case METHOD_VECTOR_CENTROID:
|
|
dist.merge(idx1, idx2, i);
|
|
break;
|
|
|
|
case METHOD_VECTOR_MEDIAN:
|
|
dist.merge_weighted(idx1, idx2, i);
|
|
break;
|
|
|
|
default:
|
|
throw std::runtime_error(std::string("Invalid method."));
|
|
}
|
|
|
|
n_nghbr[i] = active_nodes.start;
|
|
if (method==METHOD_VECTOR_WARD) {
|
|
/*
|
|
Ward linkage.
|
|
|
|
Shorter and longer distances can occur, not smaller than min(d1,d2)
|
|
but maybe bigger than max(d1,d2).
|
|
*/
|
|
min = dist.ward_extended(active_nodes.start, i);
|
|
for (j=active_nodes.succ[active_nodes.start]; j<i;
|
|
j=active_nodes.succ[j]) {
|
|
t_float tmp = dist.ward_extended(j, i);
|
|
if (tmp < min) {
|
|
min = tmp;
|
|
n_nghbr[i] = j;
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
/*
|
|
Centroid and median linkage.
|
|
|
|
Shorter and longer distances can occur, not bigger than max(d1,d2)
|
|
but maybe smaller than min(d1,d2).
|
|
*/
|
|
min = dist.sqeuclidean_extended(active_nodes.start, i);
|
|
for (j=active_nodes.succ[active_nodes.start]; j<i;
|
|
j=active_nodes.succ[j]) {
|
|
t_float tmp = dist.sqeuclidean_extended(j, i);
|
|
if (tmp < min) {
|
|
min = tmp;
|
|
n_nghbr[i] = j;
|
|
}
|
|
}
|
|
}
|
|
if (idx2<active_nodes.start) {
|
|
nn_distances.remove(active_nodes.start);
|
|
} else {
|
|
nn_distances.remove(idx2);
|
|
}
|
|
nn_distances.replace(idx1, i, min);
|
|
}
|
|
}
|
|
}
|
|
|
|
#if HAVE_VISIBILITY
|
|
#pragma GCC visibility pop
|
|
#endif
|
|
|