 /*
  * Khoros: $Id$
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id$";
#endif

 /*
  * $Log$
  */

/*
 *   Copyright, 1991, The Regents of the University of California.
 *   This software was produced under a U.S. Government contract
 *   (W-7405-ENG-36) by the Los Alamos National Laboratory, which is
 *   operated by the University of California for the U.S. Department
 *   of Energy.  The U.S. Government is licensed to use, reproduce,
 *   and distribute this software.  Neither the Government nor the
 *   University makes any warranty, express or implied, or assumes
 *   any liability responsibility for the use of this software.
 */

#include<stdio.h>

#define VALID_CLUSTER -1

/*
 *  merge_alg_1 ()    Merging Algorithm 1:  Merge clusters together based
 *                    on their Euclidean distance.  At any given iteration,
 *                    the two nearest clusters are merged.  Cluster definition
 *                    can consist of mean vector (mandatory), number of
 *                    elements, upper triangular covariance matrix, and/or
 *                    a diagonal covariance matrix.  If any of these values
 *                    are not available, the associated pointer should be 
 *                    passed in as NULL.
 *
 *                    Returned value is always 1.
 *
 *  Written by:  Patrick M. Kelly
 *  Date:        8/3/91
 */

int merge_alg_1 ( dim, num_entries, final_num, means, count, cov_upper, 
	cov_diag, mapping )

int dim;			/* Problem dimension */
int *num_entries;	/* Number of entries in the table */
int final_num;		/* Desired number of final clusters */
float **means;		/* Mean vectors */
int *count;		/* Number of elements */
float ***cov_upper;	/* Upper triangular covariance matrix */
float **cov_diag;	/* Diagonal covariance matrix */
int **mapping;		/* Mapping for the cluster reduction */

{
	int i, j, k;			/* Loop control */
	int row, col;			/* Location of smalles distance */
	float **sq_dists;		/* Squared distances */
	float max;			/* Maximum starting distance */
	int num_deleted;		/* Number of clusters deleted */
	float mindist;			/* Minimum squared distance */
	int num_cur_entries;	/* Number of current entries in set of clusters */

	float ftemp1;			/* Temporary */

	extern float **SM_get_matrix ();
	extern void UT_find_min_pos_2 ();

	/*  INITIALIZE TABLE OF SQUARED DISTANCES  */

	sq_dists = SM_get_matrix ( *num_entries );

	max = 0.0;

	for ( i = 0 ; i < *num_entries ; i ++ ) {

		sq_dists [i][i] = VALID_CLUSTER;

		for ( j = i+1 ; j < *num_entries ; j ++ ) {

			sq_dists [i][j] = 0.0;
			for ( k = 0 ; k < dim ; k ++ ) {
				ftemp1 = means [i][k] - means [j][k];
				ftemp1 *= ftemp1;
				sq_dists [i][j] += ftemp1;
			}

			if ( sq_dists [i][j] > max ) max = sq_dists [i][j] ;

		}

	}

	/*  REDUCE THE NUMBER OF CLUSTERS  */

	num_deleted = 0;

	while ( num_deleted < (*num_entries - final_num) ) {

		/*  FIND MINIMUM DISTANCE  */

		UT_find_min_pos_2 (sq_dists, *num_entries, &mindist, &row, &col, max);

		/*  MERGE CLUSTERS  */

		merge ( dim, means, count, cov_upper, cov_diag, row, col );

		/*  UPDATE DISTANCES  */

		for ( i = 0 ; i < row ; i ++ ) {
			sq_dists [i][row] = 0.0;
			for ( k = 0 ; k < dim ; k ++ ) {
				ftemp1 = means [i][k] - means [row][k];
				ftemp1 *= ftemp1;
				sq_dists [i][row] += ftemp1;
			}
		}

		for ( i = row+1 ; i < *num_entries ; i ++ ) {
			sq_dists [row][i] = 0.0;
			for ( k = 0 ; k < dim ; k ++ ) {
				ftemp1 = means [i][k] - means [row][k];
				ftemp1 *= ftemp1;
				sq_dists [row][i] += ftemp1;
			}
		}

		/*  DELETE CLUSTER BY SETTING DIAGONAL ELEMENT > 0  */

		num_deleted ++;
		sq_dists [col][col] = (float) row;    /* annexed by cluster "row" */

		/*  ANYTHING PREVIOUSLY ANNEXED BY "col" MUST BE CHANGED  */

		for ( k = 0 ; k < *num_entries ; k ++ ) {
			if ( sq_dists [k][k] == (float) col ) {
				sq_dists [k][k] = (float) row ;
			}
		}

	}

	/*  CREATE THE MAPPING FOR THE CLUSTER NUMBERS  */

	*mapping = (int *) malloc ( sizeof(int) * *num_entries );
	if ( *mapping == NULL ) {
		fprintf(stderr,"MEMORY ALLOCATION FAILURE\n");
		exit (-1);
	}

	/*  PHYSICALLY DELETE INVALID CLUSTERS  */

	num_cur_entries = 0;
	for ( i = 0 ; i < *num_entries ; i ++ ) {

		if ( sq_dists [i][i] == VALID_CLUSTER ) {

			means [num_cur_entries] = means [i];
			if ( cov_diag != NULL )
				cov_diag [num_cur_entries] = cov_diag [i];
			if ( cov_upper != NULL )
				cov_upper [num_cur_entries] = cov_upper [i];
			if ( count != NULL )
				count [num_cur_entries] = count [i];

			(*mapping) [i] = num_cur_entries ;

			num_cur_entries ++;

		} else {

			(*mapping) [i] = (*mapping) [(int)sq_dists [i][i]];

		}

		if ( (*mapping) [i] >= final_num ) 
			fprintf (stderr,"WARNING: (%d) is wrong (%d)\n",i,(*mapping)[i]);

	}

	*num_entries -= num_deleted;
	if ( *num_entries != num_cur_entries ) {
		fprintf (stderr,"WARNING:  Bogus number of entries?\n");
		fprintf (stderr,"          %d vs %d\n",*num_entries,num_cur_entries);
	}

	/*  FREE MEMORY AND EXIT  */

	free ( sq_dists [0] );
	free ( sq_dists );

     return (1);
}

