/*
 * Copyright 2009-2010 Yuichiro Moriguchi
 *
 * Licensed 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 net.morilib.awk.matrix;

import net.morilib.awk.misc.Strings;

/**
 *
 *
 * @author MORIGUCHI, Yuichiro 2011/09/19
 */
public final class DoubleMatrices {

	//
	private DoubleMatrices() {}

	/**
	 * 
	 * @param matrix
	 */
	public static void print(StringBuilder b, DoubleMatrix matrix) {
		int l = -1;
		String dlm = "";

		for(int i = 0; i < matrix.rowSize(); i++) {
			for(int j = 0; j < matrix.columnSize(); j++) {
				double d = matrix.get(i, j);
				String s = Double.toString(d);

				l = (l < s.length()) ? s.length() : l;
			}
		}

		for(int i = 0; i < matrix.rowSize(); i++) {
			b.append("[");
			dlm = "";
			for(int j = 0; j < matrix.columnSize(); j++) {
				double d = matrix.get(i, j);
				String s = Double.toString(d);

				b.append(dlm).append(Strings.rpad(s, l));
				dlm = " ";
			}
			b.append("]\n");
		}
	}

	/**
	 * 
	 * @param a
	 * @return
	 * @throws AwkMatrixException
	 */
	public static DoubleMatrix[] decomposeLU(
			DoubleMatrix a) throws AwkMatrixException {
		ArrayDoubleMatrix r = new ArrayDoubleMatrix(a);
		int[] order, order1;
		int   parity = 0;

		if(a.rowSize() != a.columnSize()) {
			throw new AwkMatrixException("square matrix required");
		}

		order = new int[a.rowSize()];
		for(int i = 0; i < a.rowSize(); i++) {
			order[i] = i;
		}

		for(int i = 0; i < a.rowSize() - 1; i++) {
			int d = i, swp;
			double mx = r.get(i, i);
			double swp2;

			for(int j = i + 1; j < a.columnSize(); j++) {
				if(mx < r.get(j, i)) {
					d = j;
					mx = r.get(j, i);
				}
			}

			if(d != i) {
				swp = order[i];
				order[i] = order[d];
				order[d] = swp;
				parity++;
				for(int j = 0; j < a.columnSize(); j++) {
					swp2 = r.get(i, j);
					r.set(i, j, r.get(d, j));
					r.set(d, j, swp2);
				}
			}

			if(r.get(i, i) == 0) {
				continue;
			} else {
				for(int j = i + 1; j < a.columnSize(); j++) {
					double x = r.get(j, i);

					x = x / r.get(i, i);
					r.set(j, i, x);
					for(int k = i + 1; k < a.rowSize(); k++) {
						double y = r.get(j, k);

						y = y - (x * r.get(i, k));
						r.set(j, k, y);
					}
				}
			}
		}

		order1 = new int[order.length];
		for(int i = 0; i < order.length; i++) {
			order1[order[i]] = i;
		}

		return new DoubleMatrix[] {
				new PermutationDoubleMatrix(order1, parity),
				new LUDecomposedLMatrix(r),
				new LUDecomposedUMatrix(r)
		};
	}

	//
	private static class Umd extends AbstractImmutableDoubleMatrix {

		//
		private DoubleMatrix a;

		//
		private Umd(DoubleMatrix a) {
			this.a = a;
		}

		public double get(int row, int column) {
			return a.get(row, column);
		}

		public int rowSize() {
			return a.rowSize();
		}

		public int columnSize() {
			return a.columnSize();
		}

		public double determinant() throws AwkMatrixException {
			return a.determinant();
		}

		public DoubleMatrix inv() throws AwkMatrixException {
			return a.inv();
		}

	}

	/**
	 * 
	 * @param x
	 * @return
	 */
	public static DoubleMatrix unmodifiable(DoubleMatrix x) {
		return (x instanceof Umd) ? x : new Umd(x);
	}

	//
	private static DoubleVector mul(DoubleMatrix m,
			DoubleVector a) throws AwkMatrixException {
		double x;
		ArrayDoubleVector v;

		if(m.columnSize() != a.size()) {
			throw new AwkMatrixException(
					"cannot multiply the given matrix");
		} else {
			v = new ArrayDoubleVector(m.columnSize());
			for(int i = 0; i < m.rowSize(); i++) {
				x = 0.0;
				for(int k = 0; k < m.columnSize(); k++) {
					x = x + m.get(i, k) * a.get(k);
				}
				v.set(i, x);
			}
			return v;
		}
	}

	/**
	 * 
	 * @param a
	 * @param b
	 * @return
	 * @throws AwkMatrixException
	 */
	public static DoubleVector solve(DoubleMatrix a,
			DoubleVector b) throws AwkMatrixException {
		DoubleMatrix[] dec;
		DoubleVector x;

		if(a.rowSize() != a.columnSize()) {
			throw new AwkMatrixException("square matrix required");
		} else if(a.rowSize() != b.size()) {
			throw new AwkMatrixException(
					"size of vector is different");
		}

		// forward
		dec = decomposeLU(a);
//		x   = dec[0].inv().mul(b);
		x   = mul(dec[0].inv(), b);
		for(int i = 1; i < a.rowSize(); i++) {
			for(int j = 0; j < i; j++) {
				x.set(i, x.get(i) - (x.get(j) * dec[1].get(i, j)));
			}
		}

		// backward
		for(int i = a.rowSize() - 1; i >= 0; i--) {
			if(dec[2].get(i, i) == 0.0) {
				throw new AwkMatrixException(
						"this matrix is not regular");
			} else {
				x.set(i, x.get(i) / dec[2].get(i, i));
				for(int j = i - 1; j >= 0; j--) {
					x.set(j, x.get(j) - x.get(i) * dec[2].get(j, i));
				}
			}
		}
		return x;
	}

	/**
	 * 
	 * @param a
	 * @param order
	 * @return
	 */
	public static double determinantAllMinor(final DoubleMatrix a,
			int order) {
		final int[][] cmb = new int[1][];
		double res = 0.0;
		int len = Math.max(a.rowSize(), a.columnSize());
		AbstractImmutableDoubleMatrix
		minor = new AbstractImmutableDoubleMatrix() {

			public double get(int row, int column) {
				return a.get(cmb[0][row], cmb[0][column]);
			}

			public int rowSize() {
				return cmb[0].length;
			}

			public int columnSize() {
				return cmb[0].length;
			}

			public double determinant(
					) throws AwkMatrixException {
				DoubleMatrix[] dc;

				if(columnSize() != rowSize()) {
					throw new AwkMatrixException(
							"square matrix required");
				} else {
					dc = DoubleMatrices.decomposeLU(this);
					// det(P) * det(U)
					return dc[0].determinant() * dc[2].determinant();
				}
			}

		};

		cmb[0] = CombinationUtils.initCombination(len, order);
		do {
			res = res + minor.determinant();
		} while((cmb[0] = CombinationUtils.nextCombination(
				len, cmb[0])) != null);
		return res;
	}

}
