/* Eigen.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 20010719
 djmw 20020402 GPL header
 djmw 20020515 Modified scaling in _drawEigenvector
*/

#include "Eigen.h"
#include "NUMmachar.h"
#include "NUMlapack.h"
#include "NUM2.h"
#include "SVD.h"

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

extern machar_Table NUMfpp;

static void Graphics_ticks (Graphics g, double min, double max, int hasNumber,
	int hasTick, int hasDottedLine, int integers)
{
	double range = max - min, scale = 1, tick = min, dtick = 1;

	if (range == 0)
	{
		return;
	}
	else if (range > 1)
	{
		while (range / scale > 10) scale *= 10;
		range /= scale;
	}
	else
	{
		while (range * scale < 10) scale /= 10;
		range *= scale;
	}

	if (range < 3) dtick = 0.5; 
	dtick *= scale;	
	tick = dtick * floor (min / dtick);
	if (tick < min) tick += dtick;
	while (tick <= max)
	{
		double num = integers ? floor (tick + 0.5) : tick;
		Graphics_markBottom (g, num, hasNumber, hasTick, hasDottedLine, NULL);
		tick += dtick;
	}
}

class_methods (Eigen, Data)
	class_method_local (Eigen, destroy)
	class_method_local (Eigen, description)
	class_method_local (Eigen, copy)
	class_method_local (Eigen, equal)
	class_method_local (Eigen, writeAscii)
	class_method_local (Eigen, readAscii)
	class_method_local (Eigen, writeBinary)
	class_method_local (Eigen, readBinary)
class_methods_end

int Eigen_init (I, long numberOfEigenvalues, long dimension)
{
	iam (Eigen);
	my numberOfEigenvalues = numberOfEigenvalues;
	my dimension = dimension;
	if (! (my eigenvalues = NUMdvector (1, numberOfEigenvalues)) ||
		! (my eigenvectors = NUMdmatrix (1, numberOfEigenvalues, 
			1, dimension))) return 0;
	return 1;
}

int Eigen_initFromSquareRoot (I, double **a, long numberOfRows, long numberOfColumns)
{
	iam (Eigen);
	SVD svd;
	long i, j, k, numberOfZeroed, numberOfEigenvalues;
	
	my dimension = numberOfColumns;
	if (! (svd = SVD_create_d (a, numberOfRows, numberOfColumns))) return 0;
	
	/*
		Make sv's that are too small zero. These values occur automatically 
		when the rank of A'A < numberOfColumns. This could occur when for
		example numberOfRows <= numberOfColumns. 
		(n points in  an n-dimensional space define maximally an n-1 
		dimensional surface for which we maximally need an n-1 dimensional
		basis.)
	*/

	numberOfZeroed = SVD_zeroSmallSingularValues (svd, 0);

	numberOfEigenvalues = numberOfColumns - numberOfZeroed;
	
	if (! Eigen_init (me, numberOfEigenvalues, numberOfColumns)) goto end;

	for (k = 0, i = 1; i <= numberOfColumns; i++)
	{
		double t = svd -> d[i];
		if (t > 0)
		{
			my eigenvalues[++k] = t * t;
			for (j = 1; j <= numberOfColumns; j++)
			{
				my eigenvectors[k][j] = svd -> v[j][i];
			}
		}
	}
	
	Eigen_sort (me);
	
end:
	forget (svd);
	return ! Melder_hasError ();
}

int Eigen_initFromSquareRootPair (I, double **a, long numberOfRows,
	long numberOfColumns, double **b, long numberOfRows_b)
{
	iam (Eigen); long i, j = 0, k, ll, n;
	double *alpha = NULL, *beta = NULL, **q = NULL, maxsv2 = -10;
	int invertR = 0, productsvd = 0;

/*	Melder_assert (numberOfRows >= numberOfColumns || numberOfRows_b >= numberOfColumns);*/
	
	my dimension = numberOfColumns;
	
	if (! (alpha = NUMdvector (1, numberOfColumns)) ||
		! (beta = NUMdvector (1, numberOfColumns)) ||
		! (q = NUMdmatrix (1, numberOfColumns, 1, numberOfColumns))) goto end;
	
	/*
		Use generalized svd to solve A'A x - lambda B'B x = 0.
	*/
	
	if (! NUMgsvdcmp (a, numberOfRows, numberOfColumns, b, numberOfRows_b,
			productsvd, &k, &ll, alpha, beta, NULL, NULL, q, invertR)) goto end;

	/*
		Don't copy eigenvalues that are < eps * maximumEigenvalue.
	*/
	
	if (! NUMfpp) NUMmachar ();
	
	for (i = 1; i <= ll; i++)
	{
		double t = alpha[k+i] / beta[k+i];
		alpha[k+i] = t * t;	
		if (alpha[k+i] > maxsv2) maxsv2 = alpha[k+i];
	}
	
	for (n = 0, i = 1; i <= ll; i++)
	{
		if (alpha[k+i] < NUMfpp -> eps_d * maxsv2)
		{
			n++; alpha[k+i] = -1;
		}
	}
		
	if (! Eigen_init (me, ll - n, numberOfColumns)) goto end;

	for (n = 0, i = 1; i <= ll; i++)
	{
		if (alpha[k+i] > 0)
		{
			n++; my eigenvalues[n] = alpha[k+i];
			for (j = 1; j <= numberOfColumns; j++)
			{
				my eigenvectors[n][j] = q[j][k+i];
			}
		}
	}
	
	NUMnormalizeRows_d (my eigenvectors, my numberOfEigenvalues, numberOfColumns, 1);
	
	Eigen_sort (me);
		 
end:
	NUMdvector_free (alpha, 1);
	NUMdvector_free (beta, 1);
	NUMdmatrix_free (q, 1, 1);
	return ! Melder_hasError();
}

int Eigen_initFromSymmetricMatrix (I, double **a, long n)
{
	iam (Eigen); double **ac = NULL, *e = NULL, *ev = NULL; long i, j;
	
	my dimension = my numberOfEigenvalues = n;
	if (! (ac = NUMdmatrix_copy (a, 1, n, 1, n)) ||
		! (ev = NUMdvector (1, n)) ||
		! (e = NUMdvector (1, n)) ||
		! Eigen_init (me, n, n)) goto end;
		
	NUMtred2_d (ac, n, ev, e);
	
	if (! NUMtqli_d (ev, e, n, ac)) goto end;
		
	for (i = 1; i <= n; i++)
	{
		my eigenvalues[i] = ev[i];
		for (j = 1; j <= n; j++) my eigenvectors[i][j] = ac[j][i];
	}
	
	Eigen_sort (me);
	
end:
	NUMdvector_free (e, 1);
	NUMdvector_free (ev, 1);
	NUMdmatrix_free (ac, 1, 1);
	return ! Melder_hasError ();
}

Eigen Eigen_create (long numberOfEigenvalues, long dimension)
{
	Eigen me = new (Eigen);
	
	if (! me || ! Eigen_init (me, numberOfEigenvalues, dimension)) forget (me);
	return me;
}

long Eigen_getNumberOfEigenvectors (I)
{
	iam (Eigen);
	return my numberOfEigenvalues;
}

long Eigen_getDimensionOfComponents (I)
{
	iam (Eigen);
	return my dimension;
}

double Eigen_getSumOfEigenvalues (I)
{
	iam (Eigen); long i; double sum = 0;
	for (i = 1; i <= my numberOfEigenvalues; i++) sum += my eigenvalues[i];
	return sum;
}

double Eigen_getCumulativeContributionOfComponents (I, long from, long to)
{
	iam (Eigen);
	long i;
	double partial = 0, sum = 0;
	
	if (from > 0 && to <= my numberOfEigenvalues && from <= to)
	{
		for (i = 1; i <= my numberOfEigenvalues; i++)
		{
			sum += my eigenvalues[i];
			if (i >= from && i <= to) partial += my eigenvalues[i];
		}
	}
	return sum > 0 ? partial / sum : 0;

}

long Eigen_getDimensionOfFraction (I, double fraction)
{
	iam (Eigen); 
	long n;
	double p, sum = Eigen_getSumOfEigenvalues (me);
	
	if (sum == 0) return 1;
	
	n = 1; p = my eigenvalues[1];
	while (p / sum < fraction && n < my numberOfEigenvalues)
	{
		p += my eigenvalues[++n];
	}
	return n;
}

void Eigen_sort (I)
{
	iam (Eigen); 
	long i, j, k;
	double *e = my eigenvalues, **v = my eigenvectors;
			
    for (i = 1; i < my numberOfEigenvalues; i++)
    {
        double emax = e[k=i];
        for (j = i + 1; j <= my numberOfEigenvalues; j++)
		{
			if (e[j] > emax) emax = e[k=j];
		}
		if (k != i)
        {
			/*
				Swap eigenvalues and eigenvectors
			*/
            e[k] = e[i]; e[i] = emax;
            for (j = 1; j <= my dimension; j++)
			{
				double tmp = v[i][j]; v[i][j] = v[k][j]; v[k][j] = tmp;
			}
        }
    }
}

void Eigen_invertEigenvector (I, long ivec)
{
	iam (Eigen); long j;

	if (ivec < 1 || ivec > my numberOfEigenvalues) return;

	for (j = 1; j <= my dimension; j++)
	{
		my eigenvectors[ivec][j] = - my eigenvectors[ivec][j];
	}
}

void Eigen_drawEigenvalues_scree (I, Graphics g, long first, long last, 
	double minimum, double maximum, double size_mm, const char *mark, 
	int garnish)
{
	iam (Eigen); double xmin = first, xmax = last; long i;

	if (last <= first )
	{
		first = 1; last = my numberOfEigenvalues;
		xmin = 0; xmax = last + 1;
	}
	if (maximum <= minimum)
	{
		maximum = my eigenvalues[first]; minimum = my eigenvalues[last];
	}
	Graphics_setInner (g);
	Graphics_setWindow (g, xmin, xmax, minimum, maximum);
	for (i=first; i <= last; i++)
	{
		Graphics_mark (g, i, my eigenvalues[i], size_mm, mark);
	}
	Graphics_unsetInner (g);
	if (garnish)
	{
    	Graphics_drawInnerBox (g);
    	Graphics_marksLeft (g, 2, 1, 1, 0);
    	Graphics_textLeft (g, 1, "eigenvalue");
		Graphics_ticks (g, first, last, 1, 1, 0, 1);
		Graphics_textBottom (g, 1, "eigenvalue number");
	}
}

void Eigen_drawEigenvector (I, Graphics g, long ivec, long first, long last,
	double minimum, double maximum, int weigh, double size_mm, const char *mark,
	int connect, int garnish)
{
	iam (Eigen); 
	long i; 
	double xmin = first, xmax = last, *vec, w;
	
	if (ivec < 1 || ivec > my numberOfEigenvalues) return;
	
	if (last <= first )
	{
		first = 1; last = my dimension;
		xmin = 0; xmax = last + 1;
	}
	vec = my eigenvectors[ivec];
	w = weigh ? sqrt (my eigenvalues[ivec]) : 1;
	if (maximum <= minimum)
	{
		NUMdvector_extrema (vec, first, last, &minimum, &maximum);
		maximum *= w; minimum *= w;
	}
	Graphics_setInner (g);
	Graphics_setWindow (g, xmin, xmax, minimum, maximum);
	
	Graphics_mark (g, first, w * vec[first], size_mm, mark);
	for (i = first + 1; i <= last; i++)
	{
		Graphics_mark (g, i, w * vec[i], size_mm, mark);
		if (connect) Graphics_line (g, i - 1, w * vec[i-1], i, w * vec[i]);
	}
	Graphics_unsetInner (g);
	if (garnish)
	{
		char ltext[20]; 
		sprintf (ltext, "eigenvector %d", ivec);
    	Graphics_drawInnerBox (g);
    	if (minimum * maximum < 0) Graphics_markLeft (g, 0.0, 1, 1, 1, NULL);
    	Graphics_marksLeft (g, 2, 1, 1, 0);
    	Graphics_textLeft (g, 1, ltext);
 		Graphics_ticks (g, first, last, 1, 1, 0, 1);
		Graphics_textBottom (g, 1, "dimension");
	}
}

/* End of file Eigen.c */
