/*
 * Bernoulli Compiler
 * Copyright (c) Cornell University
 * Department of Computer Science
 * 
 * Kamen Yotov (kamen@yotov.org)
 * 
 * $Source: C:/CVS/kyotov/kyotov/Research/BC/ILP/Math/Matrix.cs,v $
 * $Revision: 1.3 $
 * $Date: 2003/03/28 19:57:33 $
 */

using System;

namespace ILP
{
	namespace Math
	{
		public enum MatrixType
		{
			Zero,
			Identity
		}

		public class Matrix: Base2D
		{
			protected Element[,] e;
			protected Rows rows;
			protected Columns columns;

			public Matrix ()
			{
				rows = new Rows(this);
				columns = new Columns(this);
			}

			public Matrix (int r, int c) :
				this(r, c, MatrixType.Zero)
			{
			}

			public Matrix (int r, int c, MatrixType t) :
				this()
			{
				e = new Element[r, c];

				switch (t)
				{
					case MatrixType.Zero:
					{
						for (int ri = 0; ri < Rows.Size; ri++)
							for (int ci = 0; ci < Columns.Size; ci++)
								this[ri, ci] = 0;

						break;
					}
					case MatrixType.Identity:
					{
						for (int ri = 0; ri < Rows.Size; ri++)
							for (int ci = 0; ci < Columns.Size; ci++)
								this[ri, ci] = ri == ci ? 1 : 0;

						break;
					}
				}
			}

			public Matrix (int n) :
				this(n, MatrixType.Zero)
			{
			}

			public Matrix (int n, MatrixType t) :
				this(n, n, t)
			{
			}

			public Matrix (Base1D v) :
				this(new Vectors(v))
			{
			}

			public Matrix (Base2D v) :
				this()
			{
				e = new Element[v.Size, v.ElementSize];

				for (int ri = 0; ri < v.Size; ri++)
					for (int ci = 0; ci < v.ElementSize; ci++)
						this[ri, ci] = v[ri, ci];
			}

			public override int Size
			{
				get
				{
					return e != null ? e.GetUpperBound(0) + 1 : 0;
				}
			}

			public override int ElementSize
			{
				get
				{
					return e != null ? e.GetUpperBound(1) + 1 : 0;
				}
			}

			public override Base1D this [int ri]
			{
				get
				{
					return Rows[ri];
				}
				set
				{
					Rows[ri] = value;
				}
			}

			public override Element this [int ri, int ci]
			{
				get
				{
					return e[ri, ci];
				}
				set
				{
					e[ri, ci] = value;
				}
			}

			public Base2D Rows
			{
				get
				{
					return rows;
				}
				set
				{
					if (Rows.Size != value.Size)
						throw new InvalidOperationException();

					for (int ri = 0; ri < Rows.Size; ri++)
						Rows[ri] = value[ri];
				}
			}

			public Base2D Columns
			{
				get
				{
					return columns;
				}
				set
				{
					if (Columns.Size != value.Size)
						throw new InvalidOperationException();

					for (int ci = 0; ci < Columns.Size; ci++)
						Columns[ci] = value[ci];
				}
			}

			public Matrix Transposed
			{
				get
				{
					Matrix R = new Matrix(Columns.Size, Rows.Size);

					for (int i = 0; i < Rows.Size; i++)
						R.Columns[i] = Rows[i];

					return R;
				}
			}

			public Matrix Inverse
			{
				get
				{
					Matrix A = new Matrix(this);

					if (!A.IsSquare)
						throw new InvalidOperationException();

					Matrix R = new Matrix(Size, MatrixType.Identity);

					for (int i = 0; i < Size; i++)
					{
						if (A[i, i] == 0)
							for (int j = i + 1; j < Size; j++)
								if (A[j, i] != 0)
								{
									R.Rows.Swap(i, j);
									A.Rows.Swap(i, j);

									break;
								}
					
						if (A[i, i] == 0)
							throw new InvalidOperationException();

						for (int j = 0; j < Size; j++)
							if (i != j && A[j, i] != 0)
							{
								Element g = Element.gcd(A[i, i].Abs, A[j, i].Abs);

								R.Rows[j] *= A[i, i].Abs / g;
								A.Rows[j] *= A[i, i].Abs / g;

								Element m = -A[j, i] / A[i, i];

								R.Rows[j] += R.Rows[i] * m;
								A.Rows[j] += A.Rows[i] * m;
							}
					}

					for (int i = 0; i < Size; i++)
					{
						R.Rows[i] /= A[i, i];
						A.Rows[i] /= A[i, i];
					}

					return R;
				}
			}

			public Vectors Kernel
			{
				get
				{
					Matrix A = new Matrix(this);
					Matrix U = Matrix.ColumnEchelonForm(ref A);

					Vectors r = new Vectors();

					for (int ci = A.Columns.Size - 1; ci >= 0; ci--)
						if (A.Columns[ci].IsZero)
							r.Add(Vector.One(A.Columns.Size, ci));

					try
					{
						Matrix m = U * new Matrix(r).Transposed;

						return new Vectors(m.Transposed);
					}
					catch
					{
						return new Vectors();
					}
				}
			}

			public bool IsSquare
			{
				get
				{
					return Rows.Size == Columns.Size;
				}
			}

			public Element Determinant
			{
				get
				{
					if (!IsSquare)
						throw new InvalidOperationException();

					// to be implemented;

					return 1;
				}
			}

			public bool IsUnimodular
			{
				get
				{
					return Determinant.Abs == 1;
				}
			}

			public static Matrix ColumnEchelonForm (ref Matrix A)
			{
				Matrix U = new Matrix(A.Columns.Size, MatrixType.Identity);

				for (int i = 0, k = 0; i < A.Rows.Size && k < A.Columns.Size; i++)
				{
					for (int j = k + 1; j < A.Columns.Size; j++)
					{
						Element a = A[i, k];
						Element b = A[i, j];

						if (a == 0 && b == 0)
							continue;

						Element u, v, d = Element.gcd(a, b, out u, out v);

						Matrix T = new Matrix(A.Columns.Size, MatrixType.Identity);

						T[k, k] = u;
						T[j, k] = v;
						T[k, j] = -b / d;
						T[j, j] = a/d;

						A = A * T;
						U = U * T;
					}

					if (A[i, k] != 0)
						k++;
				}

				return U;
			}

			public static Matrix operator * (Matrix A, Matrix B)
			{
				if (A.Columns.Size != B.Rows.Size)
					throw new InvalidOperationException();

				Matrix r = new Matrix(A.Rows.Size, B.Columns.Size);

				for (int ri = 0; ri < r.Rows.Size; ri++)
					for (int ci = 0; ci < r.Columns.Size; ci++)
						r[ri, ci] = A.Rows[ri] ^ B.Columns[ci];

				return r;
			}

			public static Vector operator * (Matrix A, Vector v)
			{
				if (A.Columns.Size != v.Size)
					throw new InvalidOperationException();

				Vector r = new Vector(A.Rows.Size);

				for (int i = 0; i < A.Rows.Size; i++)
					r[i] = A.Rows[i] ^ v;

				return r;
			}

			public static Vectors operator * (Matrix A, Vectors v)
			{
				Vectors r = new Vectors();

				foreach (Vector t in v)
					r.Add(A * t);

				return r;
			}
		}
	}
}