 /*
  * 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.
 */


/*
 *  SM_comb_stats ()  Combine the statistics for two clusters of data.
 *
 *                   Cluster 3 <---  Cluster 1 + Cluster 2
 *
 *       (NOTE:  It is OK to call this routine with the same parameter
 *               in more than one place.  For example, you can add
 *               cluster 1 data "into" cluster 2)
 *
 *  Written by:  Patrick M. Kelly
 *  Date:        11/5/90
 *
 *  MODIFICATION HISTORY:
 *
 *       11/16/90  Patrick M. Kelly
 *                 Re-arranged factors to reduce numerical round-off
 *                 errors.
 */

void SM_comb_stats ( dim, num_1, num_2, num_3, m1, m2, m3, C1, C2, C3 )

int dim;						/* Problem Dimension */
int num_1, num_2, *num_3;		/* Number of data points */
float *m1, *m2, *m3;			/* Mean vectors */
float **C1, **C2, **C3;			/* SM covariance matrices */

{
	int i, j;					/* Loop Control */
	float num_1_f;				/* Floating point rep. of num_1 */
	float num_2_f;				/* Floating point rep. of num_2 */
	float l1, l2, l3, l4, l5;	/* Temporary */

	num_1_f = (float) num_1;
	num_2_f = (float) num_2;

	l1 = num_1_f - 1.0 ;
	l2 = num_2_f - 1.0 ;
	l5 = num_1_f + num_2_f;
	l3 = num_1_f * num_2_f / l5 ;
	l4 = num_1_f + num_2_f - 1.0 ;

	/*
	 *  Create the new covariance matrix.
	 *
	 *  (NOTE:  REMEMBER TO ONLY ACCESS THE UPPER TRIANGLE ELEMENTS
	 *          OF THE MATRICES!  SEE SM_MATRIX ROUTINES)
	 */
	for ( i = 0 ; i < dim ; i ++ ) {
	for ( j = i ; j < dim ; j ++ ) {

		C3 [i][j] = l1 * C1 [i][j] + l2 * C2 [i][j] +
				l3 * (m1 [i] - m2 [i]) * (m1 [j] - m2 [j]);

		C3 [i][j] /= l4 ;

	}
	}

	/*
	 *  Create the new mean vector.
	 */
	for ( i = 0 ; i < dim ; i ++ ) {

		m3 [i] =  ( num_1_f * m1 [i] + num_2_f  * m2 [i] ) / l5;

	}

	/*
	 *  Set number of data points in the cluster.
	 */
	*num_3 = num_1 + num_2;

	return;
}
