/* PCA.c
 * 
 * Copyright (C) 1997-2002 David Weenink
 * 
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
/* 
 djmw 20000511 PCA_and_TableOfReal_to_Configuration: added centralize option.
 	(later removed)
 djmw 20020327 PCA_and_TableOfReal_to_Configuration modified internals.
 djmw 20020418 Removed some causes for compiler warnings
 djmw 20020502 modified call Eigen_and_TableOfReal_project_into
*/

#include "PCA.h"
#include "NUMlapack.h"
#include "NUM2.h"
#include "TableOfReal_extensions.h"
#include "Eigen_and_TableOfReal.h"
#include "Configuration.h"

#include "oo_DESTROY.h"
#include "PCA_def.h"
#include "oo_COPY.h"
#include "PCA_def.h"
#include "oo_EQUAL.h"
#include "PCA_def.h"
#include "oo_WRITE_ASCII.h"
#include "PCA_def.h"
#include "oo_READ_ASCII.h"
#include "PCA_def.h"
#include "oo_WRITE_BINARY.h"
#include "PCA_def.h"
#include "oo_READ_BINARY.h"
#include "PCA_def.h"
#include "oo_DESCRIPTION.h"
#include "PCA_def.h"

static void info (I) {
	iam (PCA);
	Melder_information (
		"Number of components = %d\n"
		"Number of dimensions = %d\n"
		"Number of observations = %d",
		my numberOfEigenvalues, my dimension, my numberOfObservations);
}

class_methods (PCA, Eigen)
	class_method_local (PCA, destroy)
	class_method_local (PCA, description)
	class_method_local (PCA, copy)
	class_method_local (PCA, equal)
	class_method_local (PCA, writeAscii)
	class_method_local (PCA, readAscii)
	class_method_local (PCA, writeBinary)
	class_method_local (PCA, readBinary)
	class_method (info)
class_methods_end


PCA PCA_create (long numberOfComponents, long dimension)
{
	PCA me = new (PCA);
	if (! me || ! Eigen_init (me, numberOfComponents, dimension) ||
		! (my labels = NUMvector (sizeof (char *), 1, dimension)) ||
		! (my centroid = NUMdvector (1, dimension))) forget (me);
	return me;
}

void PCA_setNumberOfObservations (PCA me, long numberOfObservations)
{
	my numberOfObservations = numberOfObservations;
}

long PCA_getNumberOfObservations (PCA me)
{
	return my numberOfObservations;
}

void PCA_getEqualityOfEigenvalues (PCA me, long from, long to, int conservative,
	double *probability, double *chisq, long *ndf)
{
	long i, n, r; 
	double sum = 0, sumln = 0;
	
	*probability = 1; 
	*ndf = 0; 
	*chisq = 0;
	
	if (from > 0 && to == from || to > my numberOfEigenvalues) return;
	
	if (to <= from)
	{
		from = 1; 
		to = my numberOfEigenvalues;
	}
	for (i = from; i <= to; i++)
	{
		if (my eigenvalues[i] <= 0) break;
		sum += my eigenvalues[i];
		sumln += log (my eigenvalues[i]);
	}
	if (sum == 0) return;
	r = i - from;
	n = my numberOfObservations - 1;
	if (conservative) n -= from + (r * (2 * r + 1) + 2) / (6 * r);
	
	*ndf = r * (r + 1) / 2 - 1;
	*chisq = n * (r * log (sum / r) - sumln);
	*probability = NUMchiSquareQ (*chisq, *ndf);
}

PCA TableOfReal_to_PCA (I)
{
	iam (TableOfReal); 
	PCA thee; 
	double **a = NULL;
	long i, j, m = my numberOfRows, n = my numberOfColumns;
	
	if (m < 2) return Melder_errorp ("TableOfReal_to_PCA: not enough data.");

	if (m < n) Melder_warning ("TableOfReal_to_PCA: number of rows < number "
		"of columns, not enough data?");
	
	if (NUMfrobeniusnorm (m, n, my data) == 0) return Melder_errorp 
		("TableOfReal_to_PCA: some table entries must be different from zero.");
	
	if (! (thee = new (PCA)) ||
		! (a = NUMdmatrix_copy (my data, 1, m, 1, n)) ||
		! (thy centroid = NUMdvector (1, n))) goto end;
		
	for (i = 1; i <= m; i++)
	{
		for (j = 1; j <= n; j++)
		{
			thy centroid[j] += a[i][j];
		}
	}
	for (j = 1; j <= n; j++)
	{
		thy centroid[j] /= m;
	}
	for (i = 1; i <= m; i++)
	{
		for (j = 1; j <= n; j++)
		{
			a[i][j] -= thy centroid[j];
		}
	}
	
	if (! Eigen_initFromSquareRoot (thee, a, m, n) ||
		! (thy labels = NUMvector (sizeof (char *), 1, n))) goto end;
	
	NUMstrings_copyElements (my columnLabels, thy labels, 1, n);
	
	PCA_setNumberOfObservations (thee, m);
		
	/*
		The covariance matrix C = A'A / (N-1). However, we have constructed
		the eigenstructure for A'A. This has no consequences for the 
		eigenvectors, but the eigenvalues have to be divided by (N-1).
	*/
	
	for (i = 1; i <= thy numberOfEigenvalues; i++)
	{
		thy eigenvalues[i] /= (m - 1);
	}
	
end:

	NUMdmatrix_free (a, 1, 1);
	if (Melder_hasError	()) forget (thee);
	return thee;
}

Configuration PCA_and_TableOfReal_to_Configuration (PCA me, thou, 
	long numberOfDimensions)
{
	thouart (TableOfReal); 
	Configuration him; 
		
	if (numberOfDimensions == 0) numberOfDimensions = my numberOfEigenvalues;
	
	him = Configuration_create (thy numberOfRows, numberOfDimensions);
	if (him == NULL) return NULL;

	if (! Eigen_and_TableOfReal_project_into (me, thee, 1, thy numberOfColumns,
		& him, 1, numberOfDimensions) || 
		! NUMstrings_copyElements (thy rowLabels, his rowLabels, 1, 
			thy numberOfRows) ||
		! TableOfReal_setSequentialColumnLabels (him, 0,0, "pc", 1))
			forget (him);
		
	return him;
}

/* End of file PCA.c */
