Grammar Composition System

/src/miew/matrix.cs

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Text;

#if ALGLIB
using alglib;
#endif

namespace miew.Matrix
{
	using Math = System.Math;
	using String = System.String;

	[DebuggerDisplay("{_DebugString()}")]
	public class Matrix
	{
		private double[,] _values;
		private int c_rows;
		private int c_cols;

		private String _DebugString() { return String.Format("{0} rows by {1} columns", c_rows, c_cols); }

		public Matrix(int rowCount, int columnCount)
		{
			c_rows = rowCount;
			c_cols = columnCount;
			_values = new double[c_rows, c_cols];
		}

		public Matrix(double[] v)
			: this(v.Length, v.Length)
		{
			for (int i = 0; i < v.Length; i++)
				_values[i, i] = v[i];
		}

		public Matrix(double[] v, int c_take)
			: this(c_take, c_take)
		{
			c_take = Math.Min(c_take, v.Length);
			for (int i = 0; i < c_take; i++)
				_values[i, i] = v[i];
		}

		public Matrix(double[,] m)
			: this(m.GetLength(0), m.GetLength(1))
		{
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					_values[i, j] = m[i, j];
		}

		public double this[int row, int column]
		{
			get { return _values[row, column]; }
			set { _values[row, column] = value; }
		}

		public int RowCount { get { return c_rows; } }
		public int ColumnCount { get { return c_cols; } }

		public double[] Row(int i)
		{
			double[] v = new double[ColumnCount];
			for (int j = 0; j < ColumnCount; j++)
				v[j] = _values[i, j];
			return v;
		}

		public double[] Column(int j)
		{
			double[] v = new double[RowCount];
			for (int i = 0; i < RowCount; i++)
				v[i] = _values[i, j];
			return v;
		}

		public IEnumerable<double[]> Rows
		{
			get
			{
				//int x = _values.Cast<Double>().Count();	// can access all in [,] this way
				for (int i = 0; i < RowCount; i++)
					yield return Row(i);
			}
		}

		public IEnumerable<double[]> Columns
		{
			get
			{
				for (int i = 0; i < ColumnCount; i++)
					yield return Column(i);
			}
		}

		public void CleaveVertical(int left_cols, out Matrix left, out Matrix right)
		{
			if (left_cols == ColumnCount)
			{
				left = new Matrix(this);
				right = new Matrix(RowCount, 0);
			}
			else if (left_cols == 0)
			{
				left = new Matrix(RowCount, 0);
				right = new Matrix(this);
			}
			else if (left_cols > ColumnCount)
				throw new Exception();
			else
			{
				left = new double[RowCount, left_cols];
				right = new double[RowCount, ColumnCount - left_cols];
				for (int i = 0; i < RowCount; i++)
					for (int j = 0; j < ColumnCount; j++)
						if (j < left_cols)
							left[i, j] = _values[i, j];
						else
							right[i, j - left_cols] = _values[i, j];
			}
		}

		public void CleaveHorizontal(int top_rows, out Matrix top, out Matrix bottom)
		{
			if (top_rows == RowCount)
			{
				top = new Matrix(this);
				bottom = new Matrix(0, ColumnCount);
			}
			else if (top_rows == 0)
			{
				top = new Matrix(0, ColumnCount);
				bottom = new Matrix(this);
			}
			else if (top_rows > RowCount)
				throw new Exception();
			else
			{
				top = new double[top_rows, ColumnCount];
				bottom = new double[RowCount - top_rows, ColumnCount];
				for (int i = 0; i < RowCount; i++)
					for (int j = 0; j < ColumnCount; j++)
						if (i < top_rows)
							top[i, j] = _values[i, j];
						else
							bottom[i - top_rows, j] = _values[i, j];
			}
		}

		public void SetColumnValues(int j, double[] v)
		{
			if (v.Length != RowCount)
				throw new Exception();
			for (int i = 0; i < RowCount; i++)
				_values[i, j] = v[i];
		}

		public static Matrix Identity(int size)
		{
			Matrix m = new Matrix(size, size);
			for (int i = 0; i < size; i++)
				for (int j = 0; j < size; j++)
					m[i, j] = (i == j) ? 1.0 : 0.0;
			return m;
		}

		public Matrix Transpose()
		{
			Matrix m = new Matrix(c_cols, c_rows);
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					m[j, i] = this[i, j];
			return m;
		}

		public static Matrix Add(Matrix m_l, Matrix m_r)
		{
			Debug.Assert(m_l.ColumnCount == m_r.ColumnCount);
			Debug.Assert(m_l.RowCount == m_r.RowCount);

			Matrix m = new Matrix(m_l.RowCount, m_r.ColumnCount);
			for (int i = 0; i < m_l.RowCount; i++)
				for (int j = 0; j < m_l.ColumnCount; j++)
					m[i, j] = m_l[i, j] + m_r[i, j];
			return m;
		}

		public static Matrix Subtract(Matrix m_l, Matrix m_r)
		{
			Debug.Assert(m_l.ColumnCount == m_r.ColumnCount);
			Debug.Assert(m_l.RowCount == m_r.RowCount);
			Matrix m = new Matrix(m_l.RowCount, m_r.ColumnCount);
			for (int i = 0; i < m_l.RowCount; i++)
				for (int j = 0; j < m_l.ColumnCount; j++)
					m[i, j] = m_l[i, j] - m_r[i, j];
			return m;
		}

		public static Matrix Multiply(Matrix m_l, Matrix m_r)
		{
			if (m_l.ColumnCount != m_r.RowCount)
				throw new Exception();
			Matrix m = new Matrix(m_l.RowCount, m_r.ColumnCount);
			for (int i = 0; i < m.ColumnCount; i++)
			{
				for (int j = 0; j < m_l.RowCount; j++)
				{
					double value = 0.0;
					for (int k = 0; k < m_r.RowCount; k++)
						value += m_l[j, k] * m_r[k, i];
					m[j, i] = value;
				}
			}
			return m;
		}

		public static Matrix Multiply(double left, Matrix m_r)
		{
			Matrix m = new Matrix(m_r.RowCount, m_r.ColumnCount);
			for (int i = 0; i < m.RowCount; i++)
				for (int j = 0; j < m_r.ColumnCount; j++)
					m[i, j] = left * m_r[i, j];
			return m;
		}

		public static Matrix Multiply(Matrix m_l, double right)
		{
			Matrix m = new Matrix(m_l.RowCount, m_l.ColumnCount);
			for (int i = 0; i < m_l.RowCount; i++)
				for (int j = 0; j < m_l.ColumnCount; j++)
					m[i, j] = m_l[i, j] * right;
			return m;
		}

		public static Matrix Divide(Matrix m_l, double right)
		{
			Matrix m = new Matrix(m_l.RowCount, m_l.ColumnCount);
			for (int i = 0; i < m_l.RowCount; i++)
				for (int j = 0; j < m_l.ColumnCount; j++)
					m[i, j] = m_l[i, j] / right;
			return m;
		}

		public static double[] DiagonalVector(Matrix m)
		{
			Debug.Assert(m.RowCount == m.ColumnCount);

			double[] v = new double[m.RowCount];
			for (int i = 0; i < m.RowCount; i++)
				v[i] = m[i, i];
			return v;
		}

		public static double[,] ToDoubleArray(Matrix m)
		{
			double[,] result = new double[m.RowCount, m.ColumnCount];
			for (int i = 0; i < m.RowCount; i++)
				for (int j = 0; j < m.ColumnCount; j++)
					result[i, j] = m[i, j];
			return result;
		}

		public override string ToString()
		{
			StringBuilder sb = new StringBuilder();
			for (int i = 0; i < RowCount; i++)
			{
				for (int j = 0; j < ColumnCount; j++)
					sb.Append(_values[i, j].ToString("N4").PadLeft(8) + " ");
				sb.AppendLine();
			}
			return sb.ToString();
		}

		public Matrix Trim(int r0, int c0)
		{
			if (RowCount <= r0 && ColumnCount <= c0)
				return new Matrix(this);
			double[,] _new = new double[r0, c0];
			for (int i = 0; i < r0; i++)
				for (int j = 0; j < c0; j++)
					_new[i, j] = _values[i, j];
			return new Matrix(_new);
		}

		public Matrix TakeKRows(int k) { return Trim(k, ColumnCount); }
		public Matrix TakeKColumns(int k) { return Trim(RowCount, k); }
		public Matrix TakeKRowsAndColumns(int k) { return Trim(k, k); }

		public void NormalizeColumns()
		{
			double[] lengths = new double[ColumnCount];
			for (int i = 0; i < RowCount; i++)
				for (int j = 0; j < ColumnCount; j++)
					lengths[j] += Math.Pow(_values[i, j], 2);
			for (int j = 0; j < ColumnCount; j++)
				lengths[j] = Math.Sqrt(lengths[j]);

			for (int i = 0; i < RowCount; i++)
				for (int j = 0; j < ColumnCount; j++)
					_values[i, j] /= lengths[j];
		}

		public static Matrix operator +(Matrix m_l, Matrix m_r) { return Matrix.Add(m_l, m_r); }
		public static Matrix operator -(Matrix m_l, Matrix m_r) { return Matrix.Subtract(m_l, m_r); }
		public static Matrix operator *(Matrix m_l, Matrix m_r) { return Matrix.Multiply(m_l, m_r); }
		public static Matrix operator *(double left, Matrix m_r) { return Matrix.Multiply(left, m_r); }
		public static Matrix operator *(Matrix m_l, double right) { return Matrix.Multiply(m_l, right); }
		public static Matrix operator /(Matrix m_l, double right) { return Matrix.Divide(m_l, right); }
		public static implicit operator Matrix(double[,] m) { return new Matrix(m); }
		public static implicit operator double[,](Matrix m) { return ToDoubleArray(m); }

#if ALGLIB
		/// <summary>
		/// Make alglib's SVD nicer by using the above type
		/// </summary>
		public static bool SingularValueDecomposition(Matrix A, out double[] w, out Matrix V)
		{
			w = null;
			V = null;
			double[,] u = null;
			double[,] v = null;
			if (!svd.rmatrixsvd(A, A.RowCount, A.ColumnCount, 0, 1, 2, ref w, ref u, ref v))
				return false;
			V = new Matrix(v);
			return true;
		}
#endif
	};
}