src/clustering.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2008 by Matthias Ihrke   *
 *   mihrke@uni-goettingen.de   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "clustering.h"
#include "linalg.h"
#include <time.h>
#include "eeg.h"


Clusters*    kmedoids_repeat( const Array *distmat, int K, int repeat ){
  Clusters *C;
  Clusters *Cnew;
  int i;
  double ratio, tmp;
  bool ismatrix;
  matrix_CHECKSQR( ismatrix, distmat );
  if( !ismatrix) return NULL;
  int N = distmat->size[0];

  long initial_seed=(long)time(NULL);
  OptArgList *opts=optarglist("seed=long", initial_seed );

  C = cluster_init( K, N );
  ratio=DBL_MAX;
  for( i=0; i<repeat; i++ ){
     optarglist_optarg_by_key( opts, "seed" )->data_scalar=initial_seed+i*10;
     Cnew = kmedoids( distmat, K, opts );
     tmp = 1.0;
     
     tmp =  cluster_within_scatter ( distmat, Cnew ); 
     tmp /= cluster_between_scatter( distmat, Cnew ); 

     if( tmp<ratio ){
        dprintf("Run %i: %f\n", i, tmp);
        ratio=tmp; 
        cluster_copy( C, Cnew );
     }
     cluster_free( Cnew );
  }
  optarglist_free( opts );

  return C;
}

Clusters*    kmedoids(const Array *distmat, int K, OptArgList *optargs ){
  Clusters *C, *Cnew;
  int *medoids;
  int i, j, k, r, minidx, num_not_changed;
  double sum, minsum, mindist, curdist;
  bool ismatrix;
  matrix_CHECKSQR( ismatrix, distmat );
  if( !ismatrix) return NULL;
  int N = distmat->size[0];

  /* optional arguments */
  unsigned long seed=0;
  double x;
  optarg_PARSE_SCALAR( optargs, "seed", seed, unsigned long, x );

  if( K>N ){
     errprintf("K>=N: %i>=%i\n", K, N );
     return NULL;
  }

  /* memory alloc */
  medoids = (int*) malloc( K*sizeof(int) );
  C    = cluster_init(K, N);
  Cnew = cluster_init(K, N);

  /* initialization step, random init */ 
  const gsl_rng_type * T;
  gsl_rng * randgen;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  randgen = gsl_rng_alloc (T);
  if( seed==0 )
     seed=(unsigned long)time(NULL);
  gsl_rng_set (randgen, seed );

  Array *permut=array_new2( UINT, 1, K );
  for( i=0; i<K; i++ )
     array_INDEX1( permut, uint, i )=i;
  array_shuffle( permut, seed );
  for( i=0; i<N; i++ ){
     if( i<K ){ /* first each partition gets a guy */
        r = array_INDEX1( permut, uint, i );
     } else {
        r = gsl_rng_uniform_int( randgen, K );
     }
     C->clust[r][C->n[r]] = i;
     (C->n[r])++;
  }
  array_free( permut );

  dprintf("initialized cluster\n");  
#ifdef DEBUG
  cluster_print(stderr, C); 
#endif

  /*------------------ computation ----------------------------*/
  num_not_changed = 0;
  while(num_not_changed<1){ /* iterate until convergence */
     for( i=0; i<K; i++ ){ /* minimize distance to one of trials */
        minsum = DBL_MAX;
        minidx = 0;
        for( j=0; j<C->n[i]; j++ ){
          sum = 0;
          for( k=0; k<C->n[i]; k++){
             sum += mat_IDX( distmat, C->clust[i][j], C->clust[i][k] );
          }
          
          if( sum<minsum ){
             minsum = sum;
             minidx = C->clust[i][j];
          }
        }
        medoids[i] = minidx;
     }

     /* reset cluster */
     for( i=0; i<Cnew->K; i++ ){
        Cnew->n[i] = 0;
     }
     
     for( i=0; i<N; i++ ){ /* reassign clusters */
        mindist = DBL_MAX;
        minidx = 0;
        for( j=0; j<K; j++ ){ /* look for closest center */
          curdist = mat_IDX( distmat, medoids[j], i );
          if( curdist<mindist ){
             mindist=curdist;
             minidx=j;
          }
        }
        /*dprintf("Assigning trial %i to cluster %i\n", i, minidx);*/
        Cnew->clust[minidx][Cnew->n[minidx]] = i;
        Cnew->n[minidx] += 1;
     }

     if(!cluster_compare(C, Cnew))
        num_not_changed++;
     else
        cluster_copy(C, Cnew);
  }
  /*------------------ /computation ----------------------------*/

  free(medoids);
  cluster_free(Cnew);
  gsl_rng_free (randgen);

  return C;
}

int cluster_compare(const Clusters *c1, const Clusters *c2){
  int i, j, k, l, found, el;
  int c2i, prev_c2i;

  if( c1->K != c2->K ) return 1;
  /* clusters with equal number of trials? */
  for( i=0; i<c1->K; i++ ){
     found = 0;
     for( j=0; j<c2->K; j++ ){
        if(c1->n[i] == c2->n[j]){
          found = 1;
          break;
        }
     }
     if(!found)     
        return 1;
  }

  /* deep comparison */
  for( i=0; i<c1->K; i++ ){ /* each cluster from c1 */ 
     /* dprintf("i = %i\n", i); */
     prev_c2i = c1->K+1;
     c2i = 0;
     for( j=0; j<c1->n[i]; j++ ){ /* each element from c1 */
        el = c1->clust[i][j];
        /*      dprintf("j = %i\n", j);*/
        for( k=0; k<c2->K; k++ ){
          for( l=0; l<c2->n[k]; l++ ){
             if( el == c2->clust[k][l] ){
                c2i = k;
                if(prev_c2i==c1->K+1)
                  prev_c2i = k; /* first run */
             } 
          }
        }
        if( c2->n[c2i] != c1->n[i] ){
          /*          dprintf("ab: c2->n[%i]=%i, c1->n[%i]=%i\n", c2i, c2->n[c2i], i, c1->n[i]);*/
          return 1;
        }
        if(prev_c2i!=c2i){
          /*          dprintf("ab: prev_c2i=%i, c2i=%i\n", prev_c2i, c2i); */
          return 1;
        }
     }
  }
  return 0;
}

void cluster_copy(Clusters *dest, const Clusters *src){
  int i, j;

  dest->K = src->K;
  for( i=0; i<src->K; i++ ){
     dest->n[i] = src->n[i];
     for( j=0; j<src->n[i]; j++ ){
        dest->clust[i][j] = src->clust[i][j];
     }
  }
  
}

void cluster_free(Clusters *c){
  int i;
  free( c->n );
  for( i=0; i<c->K; i++)
     free(c->clust[i]);
  free(c->clust);
  free(c);
}

Clusters* cluster_init(int K, int maxN){
  Clusters *c;
  int i;

  c = (Clusters*) malloc( sizeof(Clusters) );
  c->K=K;
  c->n = (int*) malloc( K*sizeof(int) );
  c->clust = (int**) malloc( K*sizeof(int*) );
  for( i=0; i<K; i++ ){
     c->clust[i] = (int*) malloc( maxN*sizeof(int) );
     c->n[i] = 0;
  }

  return c;
}

void cluster_print(FILE *o,const Clusters *c){
  int i, j;

  fprintf(o, "Cluster with %i partitions\n", c->K);
  for( i=0; i<c->K; i++ ){
     fprintf(o, " C[%i]=\{ ", i+1);
     for( j=0; j<c->n[i]; j++ ){
        fprintf( o, "%i ", c->clust[i][j] );
     }
     fprintf( o, "}\n" );
  }
  fprintf( o, "\n" );
}

#if 0

GapStatistic* gapstat_init( GapStatistic *g, int K, int B ){
  int i;

  if( g==ALLOC_IN_FCT ){
     g = (GapStatistic*)malloc( sizeof(GapStatistic) );
  }
  g->K = K;
  g->B = B;
  g->khat = -1;
  dprintf("allocating memory\n");
  g->gapdistr = (double*)malloc( K*sizeof(double) );
  g->sk = (double*)malloc( K*sizeof(double) );
  g->Wk= (double*)malloc( K*sizeof(double) );
  g->Wkref= (double**) malloc( B*sizeof(double*) );
  for( i=0; i<B; i++ )
     g->Wkref[i] = (double*) malloc( K*sizeof(double) );
  g->progress=NULL;

  return g;
}

void gapstat_free( GapStatistic *g ){
  int i;
  free( g->Wk );
  for( i=0; i<g->B; i++ )
     free( g->Wkref[i] );
  free( g->Wkref );
  free( g->gapdistr );
  free( g->sk );
  free( g );
}

void gapstat_calculate( GapStatistic *gap, double **X, int n, int p, 
                                VectorDistanceFunction distfunction, const double** D ){
  int i, k, b;
  Clusters *C;
  double **Xr; /* reference distribution */
  double **Dr; /* distance in ref dist */
  double l_bar; /* (1/B sum_b( log(Wkref[b]) )) */

  /* init */
  Xr = (double**) malloc( n*sizeof( double* ) );
  Dr = (double**) malloc( n*sizeof( double* ) );
  for( i=0; i<n; i++ ){
     Xr[i] = (double*) malloc( p*sizeof( double ) );
     Dr[i] = (double*) malloc( n*sizeof( double ) );
  }

  if( gap->progress ){
     gap->progress( PROGRESSBAR_INIT, (gap->B+1)*gap->K );
  }

  /* calculation */
  for( k=1; k<=gap->K; k++ ){ /* data within-scatter */
     if( gap->progress )
        gap->progress( PROGRESSBAR_CONTINUE_LONG, k );
     C = kmedoids_repeat( D, n, k, 50 );
     //  cluster_print( stderr, C );
     gap->Wk[k-1] = get_within_scatter( D, n, C );
     cluster_free( C );
  }

  for( b=0; b<gap->B; b++ ){ /* monte Carlo for ref-data*/
     for( k=1; k<=gap->K; k++ ){
        if( gap->progress ){
          gap->progress( PROGRESSBAR_CONTINUE_LONG, (b+1)*gap->K+k );
        }
        
        Xr = gap_get_reference_distribution_simple( (const double **) X, n, p, Xr );
        Dr = vectordist_distmatrix( distfunction, (const double**)Xr, n, p, Dr, NULL, NULL );
        C = kmedoids_repeat( (const double**)Dr, n, k, 50 );
        //      cluster_print(stderr, C );
        gap->Wkref[b][k-1] = get_within_scatter( (const double**)Dr, n, C );
        cluster_free( C );
     }
  }
  
  for( k=1; k<=gap->K; k++ ){ /* gap distr */
     gap->gapdistr[k-1]=0;
     for( b=0; b<gap->B; b++ ){
        gap->gapdistr[k-1] += log( gap->Wkref[b][k-1] );
     }
     l_bar = ((gap->gapdistr[k-1])/(double)gap->B);
     gap->gapdistr[k-1] = l_bar - log( gap->Wk[k-1] );

     /* need std and sk */
     gap->sk[k-1] = 0;
     for( b=0; b<gap->B; b++ ){
        gap->sk[k-1] += SQR( log( gap->Wkref[b][k-1] ) - l_bar );
     }
     gap->sk[k-1] = sqrt( gap->sk[k-1]/(double)gap->B );
     gap->sk[k-1] = gap->sk[k-1] * sqrt( 1.0 + 1.0/(double)gap->B );
  }
  
  /* find best k */
  for( k=0; k<=gap->K-1; k++ ){
     if( gap->gapdistr[k] >= ((gap->gapdistr[k+1])-(gap->sk[k+1])) ){
        gap->khat = k+1;
        break;
     }
  }
  if( gap->progress )
     gap->progress( PROGRESSBAR_FINISH, 0 );

  /* free */
  for( i=0; i<n; i++ ){
     free(Xr[i]);
     free(Dr[i]);
  }
  free( Xr );
  free( Dr );
}

void gapstat_print( FILE *out, GapStatistic *g ){
  int i, j;
  double mean; 

  fprintf(stderr, "GapStatistic:\n"
             " K=%i\n"
             " B=%i\n"
             " khat=%i\n", g->K, g->B, g->khat );
  fprintf(stderr, " Gap       = [ ");
  for(i=0; i<g->K; i++ )
     fprintf(stderr, "(%i, %.2f) ", i+1, g->gapdistr[i]);
  fprintf(stderr, "]\n");

  fprintf(stderr, " sk        = [ ");
  for(i=0; i<g->K; i++ )
     fprintf(stderr, "(%i, %.4f) ", i+1, g->sk[i]);
  fprintf(stderr, "]\n");

  fprintf(stderr, " Wk        = [ ");
  for(i=0; i<g->K; i++ )
     fprintf(stderr, "(%i, %.4f) ", i+1, g->Wk[i]);
  fprintf(stderr, "]\n");

  fprintf(stderr, " <Wkref>_B = [ ");
  for(i=0; i<g->K; i++ ){
     mean = 0;
     for(j=0; j<g->B; j++){
        mean += g->Wkref[j][i];
     }
     mean /= (double)g->B;
     fprintf(stderr, "(%i, %.4f) ", i+1, mean);
  }
  fprintf(stderr, "]\n");

}

double** gap_get_reference_distribution_simple( const double **X, int n, int p, double **Xr ){
  int i, j;
  double minf, maxf;

  for( i=0; i<p; i++ ){ /* features */
     minf=DBL_MAX;
     maxf=DBL_MIN;
     for( j=0; j<n; j++ ){        /* find min/max for feature i */
        if( X[j][i]<minf )
          minf = X[j][i];
        if( X[j][i]>maxf )
          maxf = X[j][i];
     }
     
     for( j=0; j<n; j++ ){ /* uniform random values in [minf,maxf] */
        Xr[j][i] = (drand48()*maxf)+minf;
        if( isnan(Xr[j][i]) ){
          errprintf(" Xr[%i][%i] is nan\n", j, i);
        }
     }
  }

  return Xr;
}

double** gap_get_reference_distribution_svd( const double **X, int n, int p, double **Xr ){
  int i, j;
  double minf, maxf;
  gsl_matrix *gX, *gXd, *V;
  gsl_vector *work, *S;

  /* convert to gsl */
  gX = gsl_matrix_alloc( n, p );
  gXd= gsl_matrix_alloc( n, p );
  V  = gsl_matrix_alloc( n, p );
  S  = gsl_vector_alloc( p );
  work=gsl_vector_alloc( p );

  for( i=0; i<n; i++ ){
     for( j=0; j<p; j++ ){
        gsl_matrix_set( gX, i, j, X[i][j] );
     }
  }
  gsl_matrix_memcpy( gXd, gX );

  gsl_linalg_SV_decomp( gXd, V, S, work);

  gsl_blas_dgemm(CblasNoTrans, /* matrix mult */
                      CblasNoTrans,
                      1.0, gX, V, 0.0, gXd);

  for( i=0; i<p; i++ ){ /* features */
     minf=DBL_MAX;
     maxf=DBL_MIN;
     for( j=0; j<n; j++ ){        /* find min/max for feature i */
        if( gsl_matrix_get( gXd, j, i )<minf )
          minf = gsl_matrix_get( gXd, j, i );
        if( gsl_matrix_get( gXd, j, i )>maxf )
          maxf = gsl_matrix_get( gXd, j, i );
     }
     
     for( j=0; j<n; j++ ){ /* uniform random values in [minf,maxf] */
        gsl_matrix_set( gX, j, i, ((drand48()*maxf)+minf) );
     }
  }

  gsl_matrix_transpose_memcpy( gXd, V );
  gsl_blas_dgemm(CblasNoTrans,
                      CblasNoTrans,
                      1.0, gX, gXd, 0.0, V);

  for( i=0; i<n; i++ ){
     for( j=0; j<p; j++ ){
        Xr[j][i] = gsl_matrix_get( V, i, j );
     }
  }

  gsl_matrix_free( gX );
  gsl_matrix_free( gXd );
  gsl_matrix_free( V );
  gsl_vector_free( work );
  gsl_vector_free( S );

  return Xr;
}

#endif

double   cluster_within_scatter (const Array *distmat, const Clusters *c){
  double *sums;
  double W;
  int i, j, k; 
  bool ismatrix;
  matrix_CHECKSQR( ismatrix, distmat );
  if( !ismatrix) return -1;

  //  cluster_print( stderr, c);
  sums = (double*) malloc( c->K*sizeof(double) );
  W = 0;
  int numcomp=0;
  for( i=0; i<c->K; i++ ){ /* cluster loop */
     sums[i] = 0;
     for( j=0; j<(c->n[i])-1; j++ ){
        for( k=j+1; k<c->n[i]; k++){     
          sums[i] += mat_IDX( distmat, c->clust[i][j], c->clust[i][k]);
          numcomp++;
        }
     }
     //  sums[i] /= (double)2*c->n[i];
     /* dprintf("sums[i]=%f\n", sums[i]); */
     W += sums[i];
  }
  W /= (double)2*numcomp;

  free(sums);

  return W;
}

double   cluster_between_scatter(const Array *distmat, const Clusters *c){
  double W;
  int k1, k2, i, j;
  int num_comp; 
  bool ismatrix;
  matrix_CHECKSQR( ismatrix, distmat );
  if( !ismatrix) return -1;

  W = 0.0; 
  num_comp=0;
  for( k1=0; k1<c->K-1; k1++ ){
     for( k2=k1+1; k2<c->K; k2++ ){ /* compare cluster k1 and k2 */
        for( i=0; i<c->n[k1]; i++ ){
          for( j=0; j<c->n[k2]; j++ ){
             W += mat_IDX( distmat, c->clust[k1][i], c->clust[k2][j] );
             num_comp++;
          }
        }
     }
  }
  W /= (double)2*num_comp;

  return W;
}


Dendrogram*  agglomerative_clustering(const Array *distmat, LinkageFunction distfct){
  Dendrogram **nodes, *tmp; /* at the lowest level, we have N nodes */
  double min_d, cur_d=0.0;
  int min_i=0, min_j=0;
  int num_nodes;
  int i, j;
  bool ismatrix;
  matrix_CHECKSQR( ismatrix, distmat );
  if( !ismatrix) return NULL;
  int N=distmat->size[1];

  nodes = (Dendrogram**) malloc( N*sizeof( Dendrogram* ) );

  int clustidx=0;
  for( i=0; i<N; i++ ){ /* initialize terminal nodes */
     nodes[i] = dgram_init( i, NULL, NULL );
     nodes[i]->clustnum=clustidx++;
  }
  num_nodes = N;

  while( num_nodes > 1 ){
     dprintf("nmodes = %i\n", num_nodes );
     min_d = DBL_MAX;
     for( i=0; i<num_nodes-1; i++ ){
        for( j=i+1; j<num_nodes; j++ ){ 
          cur_d = distfct( distmat, nodes[i], nodes[j] );
          if( cur_d<=min_d ){
             min_d=cur_d;
             min_i=i;
             min_j=j;
          }
        }
     }
     
     /* we know min_i and min_j, now */
     tmp = dgram_init( -1, nodes[min_i], nodes[min_j] ); /* link the two nodes */
     tmp->height=min_d;
     tmp->clustnum=clustidx++;
     nodes[min_i] = tmp; /* remove min_i */
     tmp = nodes[num_nodes-1]; /* swap with last entry in current list */
     nodes[min_j] = tmp;
     num_nodes--;
     dprintf("Linking (%i,%i), height=%f\n", min_i, min_j, min_d );
  }
  tmp = nodes[0]; /* this is the root-node now */
  free( nodes );
  return tmp;
}

void _dgram_walk_matlab( const Dendrogram *t, Array *a ){
  if( t->val>=0 ){ /* leave */
     return;
  }
  mat_IDX( a, t->clustnum, 0 ) = t->left->clustnum+1;
  mat_IDX( a, t->clustnum, 1 ) = t->right->clustnum+1;
  mat_IDX( a, t->clustnum, 2 ) = t->height; 
  _dgram_walk_matlab( t->left, a );
  _dgram_walk_matlab( t->right, a );
}
Array* dgram_to_matlab( const Dendrogram *dgram ){
  int N=dgram_num_leaves( dgram );
  Array *mat=array_new2( DOUBLE, 2, N+N-1, 3 );
  _dgram_walk_matlab( dgram, mat );
  char buf[200];
  sprintf( buf, "%i-%i,:", N, N+N-2 );
  Array *rmat=array_slice( mat, buf );
  array_free(mat);
  return rmat;
}

Dendrogram* dgram_init(int val, Dendrogram *left, Dendrogram *right ){
  Dendrogram *c;
  c = (Dendrogram*) malloc( sizeof(Dendrogram) );
  c->val = val;
  c->height=-1.0;
  c->clustnum=0;
  c->left = left;
  c->right= right;
  return c;
}


void dgram_free(Dendrogram *t){
  if( t->right==NULL && t->left==NULL ){ /* base case */
     free( t );
     return;
  } else { /* recurse into sub-trees */
     if( t->right!=NULL ){
        dgram_free( t->right );
     } 
     if( t->left!=NULL ){
        dgram_free( t->left );
     }
  }
}
void _dgram_preorder_recursive( const Dendrogram *t, int *vals, int *n ){
  if( t->left==NULL && t->right==NULL ){
     if( t->val<0 ){
        errprintf("t->val<0 (%i) even though it's a terminal node\n", t->val);
        return;
     }
     vals[(*n)++] = t->val;
     return;
  } 
  if( t->left!=NULL ){
     //  dprintf("walking  left subtree\n");
     _dgram_preorder_recursive( t->left, vals, n );
  }
  if( t->right!=NULL ){
     _dgram_preorder_recursive( t->right, vals, n );
  }
  return;
}
void         dgram_preorder( const Dendrogram *t, int *vals, int *n ){
  int m;

  m = 0;
  //  dprintf(" starting recursive walk, m=%i\n", m);
  //  dgram_print_node(t);
  _dgram_preorder_recursive( t, vals, &m );

  /* fprintf(stderr, "vals=["); */
  /* for( i=0; i<m; i++ ){ */
  /*     fprintf(stderr, "%i ", vals[i]); */
  /* } */
  /* fprintf(stderr, "]\n"); */
  //  dprintf(" done with recursive walk, m=%i\n", m);
  *n = m;

  return;
}
double dgram_dist_singlelinkage  (const Array *d, const Dendrogram *c1, const Dendrogram *c2){
  double dist;
  int *el1, *el2;
  int n1, n2;
  int i, j;
  bool ismatrix;
  matrix_CHECK( ismatrix,d );
  if( !ismatrix ) return -1;
  if( d->size[0]!=d->size[1] ){  
     errprintf("distmat must be N x N, got (%i,%i)\n",d->size[0],d->size[1]);
     return -1; 
  }
  int N = d->size[0];

  dist = DBL_MAX;
  
  el1 = (int*) malloc( N*sizeof( int ) ); /* N suffices */
  el2 = (int*) malloc( N*sizeof( int ) ); 

  dgram_preorder( c1, el1, &n1 );
  dgram_preorder( c2, el2, &n2 );

  for( i=0; i<n1; i++ ){
     for( j=0; j<n2; j++ ){
        if( el1[i]<0 || el1[i]>=N || el2[j]<0 || el2[j]>=N ){
          errprintf("something's wrong, el[i] not in range 0-%i\n", N);
          return -1;
        }
        if( mat_IDX( d, el1[i], el2[j]) < dist ){
          dist = mat_IDX( d, el1[i], el2[j]);
        }
     }
  }
  
  //dprintf("SL-distance = %f\n", dist);

  free(el1);
  free(el2);
  return dist;
}

double dgram_dist_completelinkage(const Array *d, const Dendrogram *c1, const Dendrogram *c2){
  double dist;
  int *el1, *el2;
  int n1, n2;
  int i, j;
  bool ismatrix;
  matrix_CHECK( ismatrix,d );
  if( !ismatrix ) return -1; 
  if( d->size[0]!=d->size[1] ){  
     errprintf("distmat must be N x N, got (%i,%i)\n",d->size[0],d->size[1]);
     return -1; 
  }
  int N = d->size[0];

  dist = DBL_MIN;
  
  el1 = (int*) malloc( N*sizeof( int ) ); /* N suffices */
  el2 = (int*) malloc( N*sizeof( int ) ); 

  dgram_preorder( c1, el1, &n1 );
  dgram_preorder( c2, el2, &n2 );

  for( i=0; i<n1; i++ ){
     for( j=0; j<n2; j++ ){
        if( el1[i]<0 || el1[i]>=N || el2[j]<0 || el2[j]>=N ){
          errprintf("something's wrong, el[i] not in range 0-%i\n", N);
          return -1;
        }
        if( mat_IDX( d, el1[i], el2[j] )>dist ){
          dist = mat_IDX( d, el1[i], el2[j]);
        }
     }
  }
  
  //  dprintf("SL-distance = %f\n", dist);

  free(el1);
  free(el2);
  return dist;
}


#define END_NODE 0
#define INTERMEDIATE_NODE 1
#define NULL_NODE 2
int _dgram_get_deepest_recursive( Dendrogram *c, Dendrogram **candidate, int *depth, int curdepth ){
  int status_left=NULL_NODE;
  int status_right=NULL_NODE;

  if( c->left==NULL && c->right==NULL ){
     return END_NODE;
  }

  if( c->left!=NULL ){
     status_left  = _dgram_get_deepest_recursive( c->left,  candidate, depth, curdepth+1 );  
  }
  if( c->right!=NULL ){
     status_right = _dgram_get_deepest_recursive( c->right, candidate, depth, curdepth+1 );  
  }

  if( status_left==END_NODE && status_right==END_NODE ){

     if( curdepth>=*depth ){
        *candidate = c;
        *depth = curdepth;
     }
     dprintf("END_NODE visited, cd=%i, d=%i, c=%p, candidate=%p\n", curdepth, *depth, c, candidate);
  }

  return INTERMEDIATE_NODE;
}
Dendrogram* dgram_get_deepest( Dendrogram *c ){
  Dendrogram **d, *tmp;
  int depth;
  d = (Dendrogram **) malloc( sizeof( Dendrogram* ) );

  depth = 0;
  _dgram_get_deepest_recursive( c, d, &depth, 0 );
  dprintf("deepest=%p, depth=%i\n", *d, depth);
  tmp = *d;
  free(d);
  return tmp;
}


void _dgram_print_recursive( Dendrogram *t, int level ){
  int i;

  fprintf( stdout, "'%p(%4i,%.2f,%4i)' - ", t, t->clustnum, t->height, t->val );
  if( t->left ){ /* there is a left tree, print it */
     _dgram_print_recursive( t->left, level+1 );
  }  
  /* we need to add spaces */
  fprintf( stdout, "\n" );
  for( i=0; i<level; i++)
     fprintf( stdout, "           " );

  if( t->right ){
     _dgram_print_recursive( t->right, level+1 );
  }
}
void         dgram_print( Dendrogram *t ){
  _dgram_print_recursive( t, 0 );
  fprintf( stdout, "\n");
}

int dgram_num_leaves( const Dendrogram *t ){
  int r=0;
  if( t==NULL ){
     return 0;
  } else if( t->val!=-1 ){
     r=1;
  }
  r += dgram_num_leaves( t->left );
  r += dgram_num_leaves( t->right );
  return r;
}

void dgram_print_node( Dendrogram *t ){
  fprintf( stdout, "t=%p, val=%i, height=%f, left=%p, right=%p\n", t, t->val, t->height, t->left, t->right );
}
#undef END_NODE 
#undef INTERMEDIATE_NODE 
#undef NULL_NODE

#if 0

int      eeg_best_num_clusters_gapstat( const EEG *eeg, VectorDistanceFunction distfunction, OptArgList *optargs ){
  GapStatistic *gap;
  double x;
  int max_num_clusters;
  int num_ref_dist;
  EEG *meaneeg;
  double **D;
  int bestclust=0;

  max_num_clusters = 10;
  num_ref_dist = eeg->ntrials;
  if( optarglist_has_key( optargs, "max_num_clusters" ) ){
     x = optarglist_scalar_by_key( optargs, "max_num_clusters" );
     if( !isnan( x ) ) max_num_clusters=(int)x;
  }
  if( optarglist_has_key( optargs, "num_ref_dist" ) ){
     x = optarglist_scalar_by_key( optargs, "num_ref_dist" );
     if( !isnan( x ) ) num_ref_dist=(int)x;
  }

  meaneeg = eeg_average_channels( eeg );
  eeg_print( stderr, meaneeg, 3 );
  D = eeg_distmatrix( meaneeg, distfunction, ALLOC_IN_FCT, optargs );

  gap = gapstat_init( NULL, max_num_clusters, num_ref_dist );
  gapstat_calculate( gap, meaneeg->data[0], meaneeg->ntrials, meaneeg->n, 
                            distfunction, (const double**)D );

  bestclust = gap->khat;
#ifdef DEBUG
  gapstat_print( stderr, gap );
#endif 

  dblpp_free( D, meaneeg->ntrials );
  eeg_free( meaneeg );
  gapstat_free( gap );

  return bestclust;
}
#endif