/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math.linear;

import java.util.Arrays;
import java.util.Random;

import org.apache.commons.math.linear.EigenDecomposition;
import org.apache.commons.math.linear.EigenDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.TriDiagonalTransformer;
import org.apache.commons.math.util.MathUtils;

import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;

public class EigenDecompositionImplTest extends TestCase {

	private double[] refValues;
	private RealMatrix matrix;

	public EigenDecompositionImplTest(String name) {
		super(name);
	}

	public static Test suite() {
		TestSuite suite = new TestSuite(EigenDecompositionImplTest.class);
		suite.setName("EigenDecompositionImpl Tests");
		return suite;
	}

	public void testDimension1() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
	}

	public void testDimension2() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(new double[][] { { 59.0, 12.0 }, { 12.0, 66.0 } });
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
		assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
	}

	public void testDimension3() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(new double[][] { { 39632.0, -4824.0, -16560.0 },
				{ -4824.0, 8693.0, 7920.0 }, { -16560.0, 7920.0, 17300.0 } });
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
		assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
		assertEquals(3125.0, ed.getRealEigenvalue(2), 3.0e-11);
	}

	public void testDimension4WithSplit() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(new double[][] { { 0.784, -0.288, 0.000, 0.000 },
				{ -0.288, 0.616, 0.000, 0.000 }, { 0.000, 0.000, 0.164, -0.048 }, { 0.000, 0.000, -0.048, 0.136 } });
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
		assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
		assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
		assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
	}

	public void testDimension4WithoutSplit() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(
				new double[][] { { 0.5608, -0.2016, 0.1152, -0.2976 }, { -0.2016, 0.4432, -0.2304, 0.1152 },
						{ 0.1152, -0.2304, 0.3088, -0.1344 }, { -0.2976, 0.1152, -0.1344, 0.3872 } });
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
		assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
		assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
		assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
	}

	/** test a matrix already in tridiagonal form. */
	public void testTridiagonal() {
		Random r = new Random(4366663527842l);
		double[] ref = new double[30];
		for (int i = 0; i < ref.length; ++i) {
			if (i < 5) {
				ref[i] = 2 * r.nextDouble() - 1;
			} else {
				ref[i] = 0.0001 * r.nextDouble() + 6;
			}
		}
		Arrays.sort(ref);
		TriDiagonalTransformer t = new TriDiagonalTransformer(createTestMatrix(r, ref));
		EigenDecomposition ed = new EigenDecompositionImpl(t.getMainDiagonalRef(), t.getSecondaryDiagonalRef(),
				MathUtils.SAFE_MIN);
		double[] eigenValues = ed.getRealEigenvalues();
		assertEquals(ref.length, eigenValues.length);
		for (int i = 0; i < ref.length; ++i) {
			assertEquals(ref[ref.length - i - 1], eigenValues[i], 2.0e-14);
		}

	}

	/** test dimensions */
	public void testDimensions() {
		final int m = matrix.getRowDimension();
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		assertEquals(m, ed.getV().getRowDimension());
		assertEquals(m, ed.getV().getColumnDimension());
		assertEquals(m, ed.getD().getColumnDimension());
		assertEquals(m, ed.getD().getColumnDimension());
		assertEquals(m, ed.getVT().getRowDimension());
		assertEquals(m, ed.getVT().getColumnDimension());
	}

	/** test eigenvalues */
	public void testEigenvalues() {
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		double[] eigenValues = ed.getRealEigenvalues();
		assertEquals(refValues.length, eigenValues.length);
		for (int i = 0; i < refValues.length; ++i) {
			assertEquals(refValues[i], eigenValues[i], 3.0e-15);
		}
	}

	/** test eigenvalues for a big matrix. */
	public void testBigMatrix() {
		Random r = new Random(17748333525117l);
		double[] bigValues = new double[200];
		for (int i = 0; i < bigValues.length; ++i) {
			bigValues[i] = 2 * r.nextDouble() - 1;
		}
		Arrays.sort(bigValues);
		EigenDecomposition ed = new EigenDecompositionImpl(createTestMatrix(r, bigValues), MathUtils.SAFE_MIN);
		double[] eigenValues = ed.getRealEigenvalues();
		assertEquals(bigValues.length, eigenValues.length);
		for (int i = 0; i < bigValues.length; ++i) {
			assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14);
		}
	}

	/** test eigenvectors */
	public void testEigenvectors() {
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		for (int i = 0; i < matrix.getRowDimension(); ++i) {
			double lambda = ed.getRealEigenvalue(i);
			RealVector v = ed.getEigenvector(i);
			RealVector mV = matrix.operate(v);
			assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13);
		}
	}

	/** test A = VDVt */
	public void testAEqualVDVt() {
		EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
		RealMatrix v = ed.getV();
		RealMatrix d = ed.getD();
		RealMatrix vT = ed.getVT();
		double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm();
		assertEquals(0, norm, 6.0e-13);
	}

	/** test that V is orthogonal */
	public void testVOrthogonal() {
		RealMatrix v = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN).getV();
		RealMatrix vTv = v.transpose().multiply(v);
		RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
		assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
	}

	/** test diagonal matrix */
	public void testDiagonal() {
		double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
		RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
		EigenDecomposition ed = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN);
		assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
		assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
		assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
		assertEquals(diagonal[3], ed.getRealEigenvalue(0), 2.0e-15);
	}

	/**
	 * Matrix with eigenvalues {8, -1, -1}
	 */
	public void testRepeatedEigenvalue() {
		RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] { { 3, 2, 4 }, { 2, 0, 2 }, { 4, 2, 3 } });
		EigenDecomposition ed = new EigenDecompositionImpl(repeated, MathUtils.SAFE_MIN);
		checkEigenValues((new double[] { 8, -1, -1 }), ed, 1E-12);
		checkEigenVector((new double[] { 2, 1, 2 }), ed, 1E-12);
	}

	/**
	 * Matrix with eigenvalues {2, 0, 12}
	 */
	public void testDistinctEigenvalues() {
		RealMatrix distinct = MatrixUtils
				.createRealMatrix(new double[][] { { 3, 1, -4 }, { 1, 3, -4 }, { -4, -4, 8 } });
		EigenDecomposition ed = new EigenDecompositionImpl(distinct, MathUtils.SAFE_MIN);
		checkEigenValues((new double[] { 2, 0, 12 }), ed, 1E-12);
		checkEigenVector((new double[] { 1, -1, 0 }), ed, 1E-12);
		checkEigenVector((new double[] { 1, 1, 1 }), ed, 1E-12);
		checkEigenVector((new double[] { -1, -1, 2 }), ed, 1E-12);
	}

	/**
	 * Verifies that the given EigenDecomposition has eigenvalues equivalent to
	 * the targetValues, ignoring the order of the values and allowing values to
	 * differ by tolerance.
	 */
	protected void checkEigenValues(double[] targetValues, EigenDecomposition ed, double tolerance) {
		double[] observed = ed.getRealEigenvalues();
		for (int i = 0; i < observed.length; i++) {
			assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
			assertTrue(isIncludedValue(targetValues[i], observed, tolerance));
		}
	}

	/**
	 * Returns true iff there is an entry within tolerance of value in
	 * searchArray.
	 */
	private boolean isIncludedValue(double value, double[] searchArray, double tolerance) {
		boolean found = false;
		int i = 0;
		while (!found && i < searchArray.length) {
			if (Math.abs(value - searchArray[i]) < tolerance) {
				found = true;
			}
			i++;
		}
		return found;
	}

	/**
	 * Returns true iff eigenVector is a scalar multiple of one of the columns
	 * of ed.getV(). Does not try linear combinations - i.e., should only be
	 * used to find vectors in one-dimensional eigenspaces.
	 */
	protected void checkEigenVector(double[] eigenVector, EigenDecomposition ed, double tolerance) {
		assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
	}

	/**
	 * Returns true iff there is a column that is a scalar multiple of column in
	 * searchMatrix (modulo tolerance)
	 */
	private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix, double tolerance) {
		boolean found = false;
		int i = 0;
		while (!found && i < searchMatrix.getColumnDimension()) {
			double multiplier = 1.0;
			boolean matching = true;
			int j = 0;
			while (matching && j < searchMatrix.getRowDimension()) {
				double colEntry = searchMatrix.getEntry(j, i);
				// Use the first entry where both are non-zero as scalar
				if (Math.abs(multiplier - 1.0) <= Math.ulp(1.0) && Math.abs(colEntry) > 1E-14
						&& Math.abs(column[j]) > 1e-14) {
					multiplier = colEntry / column[j];
				}
				if (Math.abs(column[j] * multiplier - colEntry) > tolerance) {
					matching = false;
				}
				j++;
			}
			found = matching;
			i++;
		}
		return found;
	}

	@Override
	public void setUp() {
		refValues = new double[] { 2.003, 2.002, 2.001, 1.001, 1.000, 0.001 };
		matrix = createTestMatrix(new Random(35992629946426l), refValues);
	}

	@Override
	public void tearDown() {
		refValues = null;
		matrix = null;
	}

	static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
		final int n = eigenValues.length;
		final RealMatrix v = createOrthogonalMatrix(r, n);
		final RealMatrix d = createDiagonalMatrix(eigenValues, n, n);
		return v.multiply(d).multiply(v.transpose());
	}

	public static RealMatrix createOrthogonalMatrix(final Random r, final int size) {

		final double[][] data = new double[size][size];

		for (int i = 0; i < size; ++i) {
			final double[] dataI = data[i];
			double norm2 = 0;
			do {

				// generate randomly row I
				for (int j = 0; j < size; ++j) {
					dataI[j] = 2 * r.nextDouble() - 1;
				}

				// project the row in the subspace orthogonal to previous rows
				for (int k = 0; k < i; ++k) {
					final double[] dataK = data[k];
					double dotProduct = 0;
					for (int j = 0; j < size; ++j) {
						dotProduct += dataI[j] * dataK[j];
					}
					for (int j = 0; j < size; ++j) {
						dataI[j] -= dotProduct * dataK[j];
					}
				}

				// normalize the row
				norm2 = 0;
				for (final double dataIJ : dataI) {
					norm2 += dataIJ * dataIJ;
				}
				final double inv = 1.0 / Math.sqrt(norm2);
				for (int j = 0; j < size; ++j) {
					dataI[j] *= inv;
				}

			} while (norm2 * size < 0.01);
		}

		return MatrixUtils.createRealMatrix(data);

	}

	public static RealMatrix createDiagonalMatrix(final double[] diagonal, final int rows, final int columns) {
		final double[][] dData = new double[rows][columns];
		for (int i = 0; i < Math.min(rows, columns); ++i) {
			dData[i][i] = diagonal[i];
		}
		return MatrixUtils.createRealMatrix(dData);
	}

}
