/*
 * 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 org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecomposition;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;

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

public class LUDecompositionImplTest extends TestCase {
	private double[][] testData = { { 1.0, 2.0, 3.0 }, { 2.0, 5.0, 3.0 }, { 1.0, 0.0, 8.0 } };
	private double[][] testDataMinus = { { -1.0, -2.0, -3.0 }, { -2.0, -5.0, -3.0 }, { -1.0, 0.0, -8.0 } };
	private double[][] luData = { { 2.0, 3.0, 3.0 }, { 0.0, 5.0, 7.0 }, { 6.0, 9.0, 8.0 } };

	// singular matrices
	private double[][] singular = { { 2.0, 3.0 }, { 2.0, 3.0 } };
	private double[][] bigSingular = { { 1.0, 2.0, 3.0, 4.0 }, { 2.0, 5.0, 3.0, 4.0 }, { 7.0, 3.0, 256.0, 1930.0 },
			{ 3.0, 7.0, 6.0, 8.0 } }; // 4th row = 1st + 2nd

	private static final double entryTolerance = 10e-16;

	private static final double normTolerance = 10e-14;

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

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

	/** test dimensions */
	public void testDimensions() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
		LUDecomposition LU = new LUDecompositionImpl(matrix);
		assertEquals(testData.length, LU.getL().getRowDimension());
		assertEquals(testData.length, LU.getL().getColumnDimension());
		assertEquals(testData.length, LU.getU().getRowDimension());
		assertEquals(testData.length, LU.getU().getColumnDimension());
		assertEquals(testData.length, LU.getP().getRowDimension());
		assertEquals(testData.length, LU.getP().getColumnDimension());

	}

	/** test non-square matrix */
	public void testNonSquare() {
		try {
			new LUDecompositionImpl(MatrixUtils.createRealMatrix(new double[3][2]));
		} catch (InvalidMatrixException ime) {
			// expected behavior
		} catch (Exception e) {
			fail("wrong exception caught");
		}
	}

	/** test PA = LU */
	public void testPAEqualLU() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
		LUDecomposition lu = new LUDecompositionImpl(matrix);
		RealMatrix l = lu.getL();
		RealMatrix u = lu.getU();
		RealMatrix p = lu.getP();
		double norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
		assertEquals(0, norm, normTolerance);

		matrix = MatrixUtils.createRealMatrix(testDataMinus);
		lu = new LUDecompositionImpl(matrix);
		l = lu.getL();
		u = lu.getU();
		p = lu.getP();
		norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
		assertEquals(0, norm, normTolerance);

		matrix = MatrixUtils.createRealIdentityMatrix(17);
		lu = new LUDecompositionImpl(matrix);
		l = lu.getL();
		u = lu.getU();
		p = lu.getP();
		norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
		assertEquals(0, norm, normTolerance);

		matrix = MatrixUtils.createRealMatrix(singular);
		lu = new LUDecompositionImpl(matrix);
		assertFalse(lu.getSolver().isNonSingular());
		assertNull(lu.getL());
		assertNull(lu.getU());
		assertNull(lu.getP());

		matrix = MatrixUtils.createRealMatrix(bigSingular);
		lu = new LUDecompositionImpl(matrix);
		assertFalse(lu.getSolver().isNonSingular());
		assertNull(lu.getL());
		assertNull(lu.getU());
		assertNull(lu.getP());

	}

	/** test that L is lower triangular with unit diagonal */
	public void testLLowerTriangular() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
		RealMatrix l = new LUDecompositionImpl(matrix).getL();
		for (int i = 0; i < l.getRowDimension(); i++) {
			assertEquals(l.getEntry(i, i), 1, entryTolerance);
			for (int j = i + 1; j < l.getColumnDimension(); j++) {
				assertEquals(l.getEntry(i, j), 0, entryTolerance);
			}
		}
	}

	/** test that U is upper triangular */
	public void testUUpperTriangular() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
		RealMatrix u = new LUDecompositionImpl(matrix).getU();
		for (int i = 0; i < u.getRowDimension(); i++) {
			for (int j = 0; j < i; j++) {
				assertEquals(u.getEntry(i, j), 0, entryTolerance);
			}
		}
	}

	/** test that P is a permutation matrix */
	public void testPPermutation() {
		RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
		RealMatrix p = new LUDecompositionImpl(matrix).getP();

		RealMatrix ppT = p.multiply(p.transpose());
		RealMatrix id = MatrixUtils.createRealIdentityMatrix(p.getRowDimension());
		assertEquals(0, ppT.subtract(id).getNorm(), normTolerance);

		for (int i = 0; i < p.getRowDimension(); i++) {
			int zeroCount = 0;
			int oneCount = 0;
			int otherCount = 0;
			for (int j = 0; j < p.getColumnDimension(); j++) {
				final double e = p.getEntry(i, j);
				if (e == 0) {
					++zeroCount;
				} else if (e == 1) {
					++oneCount;
				} else {
					++otherCount;
				}
			}
			assertEquals(p.getColumnDimension() - 1, zeroCount);
			assertEquals(1, oneCount);
			assertEquals(0, otherCount);
		}

		for (int j = 0; j < p.getColumnDimension(); j++) {
			int zeroCount = 0;
			int oneCount = 0;
			int otherCount = 0;
			for (int i = 0; i < p.getRowDimension(); i++) {
				final double e = p.getEntry(i, j);
				if (e == 0) {
					++zeroCount;
				} else if (e == 1) {
					++oneCount;
				} else {
					++otherCount;
				}
			}
			assertEquals(p.getRowDimension() - 1, zeroCount);
			assertEquals(1, oneCount);
			assertEquals(0, otherCount);
		}

	}

	/** test singular */
	public void testSingular() {
		LUDecomposition lu = new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
		assertTrue(lu.getSolver().isNonSingular());
		lu = new LUDecompositionImpl(MatrixUtils.createRealMatrix(singular));
		assertFalse(lu.getSolver().isNonSingular());
		lu = new LUDecompositionImpl(MatrixUtils.createRealMatrix(bigSingular));
		assertFalse(lu.getSolver().isNonSingular());
	}

	/** test matrices values */
	public void testMatricesValues1() {
		LUDecomposition lu = new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
		RealMatrix lRef = MatrixUtils
				.createRealMatrix(new double[][] { { 1.0, 0.0, 0.0 }, { 0.5, 1.0, 0.0 }, { 0.5, 0.2, 1.0 } });
		RealMatrix uRef = MatrixUtils
				.createRealMatrix(new double[][] { { 2.0, 5.0, 3.0 }, { 0.0, -2.5, 6.5 }, { 0.0, 0.0, 0.2 } });
		RealMatrix pRef = MatrixUtils
				.createRealMatrix(new double[][] { { 0.0, 1.0, 0.0 }, { 0.0, 0.0, 1.0 }, { 1.0, 0.0, 0.0 } });
		int[] pivotRef = { 1, 2, 0 };

		// check values against known references
		RealMatrix l = lu.getL();
		assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
		RealMatrix u = lu.getU();
		assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
		RealMatrix p = lu.getP();
		assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
		int[] pivot = lu.getPivot();
		for (int i = 0; i < pivotRef.length; ++i) {
			assertEquals(pivotRef[i], pivot[i]);
		}

		// check the same cached instance is returned the second time
		assertTrue(l == lu.getL());
		assertTrue(u == lu.getU());
		assertTrue(p == lu.getP());

	}

	/** test matrices values */
	public void testMatricesValues2() {
		LUDecomposition lu = new LUDecompositionImpl(MatrixUtils.createRealMatrix(luData));
		RealMatrix lRef = MatrixUtils
				.createRealMatrix(new double[][] { { 1.0, 0.0, 0.0 }, { 0.0, 1.0, 0.0 }, { 1.0 / 3.0, 0.0, 1.0 } });
		RealMatrix uRef = MatrixUtils
				.createRealMatrix(new double[][] { { 6.0, 9.0, 8.0 }, { 0.0, 5.0, 7.0 }, { 0.0, 0.0, 1.0 / 3.0 } });
		RealMatrix pRef = MatrixUtils
				.createRealMatrix(new double[][] { { 0.0, 0.0, 1.0 }, { 0.0, 1.0, 0.0 }, { 1.0, 0.0, 0.0 } });
		int[] pivotRef = { 2, 1, 0 };

		// check values against known references
		RealMatrix l = lu.getL();
		assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
		RealMatrix u = lu.getU();
		assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
		RealMatrix p = lu.getP();
		assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
		int[] pivot = lu.getPivot();
		for (int i = 0; i < pivotRef.length; ++i) {
			assertEquals(pivotRef[i], pivot[i]);
		}

		// check the same cached instance is returned the second time
		assertTrue(l == lu.getL());
		assertTrue(u == lu.getU());
		assertTrue(p == lu.getP());

	}

}
