Matrix.kt

import java.util.*

/**
 * Represents a rectangular matrix of doubles
 * Created by Sergey on 19.02.2015.
 */
data class Matrix(val rows: List<List<Double>>) {

    constructor(vararg rows: List<Double>): this(rows.toList())

    init {
        if (rows.size == 0)
            throw IllegalArgumentException("Matrix components should have non-zero length.")
        if (rows.first().size.let { sz -> rows.any { it.size != sz } })
            throw IllegalArgumentException("All of the matrix rows should have the same length.")
    }

    /**
     * Checks if all the components of this and another matrix differ not more than by epsilon.
     * @param that    A matrix to compare with.
     * @param epsilon Permitted difference.
     * @return For all i, j: |this[i, j] - that[i, j]| <= epsilon
     */
    fun similar(that: Matrix, epsilon: Double = SIMILARITY_EPSILON): Boolean {
        if (height != that.height || width != that.width) return false
        return rowIndices.all { i ->
            colIndices.all { j ->
                Math.abs(this[i, j] - that[i, j]) <= epsilon
            }
        }
    }

    /**
     * Calculates the matrix m-norm.
     * @return ||A|| = max {1..n} { Sum {1..m} |a_ij| }
     */
    fun norm(): Double {
        return rows.map { it.max()!! }.max()!!
    }

    /**
     * Calculates the matrix conditionality value.
     * @return Cond A = ||A|| * ||A^-1||
     */
    fun conditionality(): Double {
        return norm() * reversed()!!.norm()
    }

    /**
     * Returns an item at i-th row, j-th column
     * @param i Row
     * @param j Column
     * @return Matrix item at i-th row, j-th column
     */
    operator fun get(i: Int, j: Int): Double {
        return rows[i][j]
    }

    val width: Int get() = rows[0].size

    val height: Int get() = rows.size

    public val rowIndices: IntRange get() = rows.indices

    public val colIndices: IntRange get() = rows[0].indices

    /**
     * Adds this matrix to another one.
     * @param that Matrix to add this to.
     * @return A = this + m
     */
    operator fun plus(that: Matrix): Matrix {
        if (height != that.height || width != that.width)
            throw IllegalArgumentException("Matrix should have the same size to be added.")
        return Matrix(
                rowIndices.map { i ->
                    colIndices.map { j -> this[i, j] + that[i, j] }
                })
    }

    operator fun minus(b: Matrix): Matrix {
        return this + b * -1.0
    }

    /**
     * Multiplies this matrix by another one.
     * @param m Matrix to multiply this by.
     * @return A = this * m;
     */
    operator fun times(m: Matrix): Matrix {
        if (width != m.height)
            throw IllegalArgumentException("Matrices cannot be multiplied due to their dimensions.")
        return Matrix(
                rowIndices.map { i ->
                    m.colIndices.map { j ->
                        colIndices.map { k -> get(i, k) * m[k, j] }.sum()
                    }
                })
    }

    /**
     * Multiplies this matrix by number.
     * @param a Number to multiply this matrix by.
     * @return A = this * a
     */
    operator fun times(a: Double): Matrix = Matrix(
            rowIndices.map { i ->
                colIndices.map { j -> this[i, j] * a }
            })

    /**
     * Multiplies the matrix by vector v as Av.
     * @param v Vector to multiply the matrix by.
     * @return Vector y = Av.
     */
    operator fun times(v: Vector): Vector {
        if (width != v.dimensions)
            throw IllegalArgumentException("Vector cannot be multiplied due to its dimensions.")
        return Vector(rowIndices.map { i ->
            colIndices.map { j -> this[i, j] * v[j] }.sum()
        })
    }

    /**
     * Get an array of the matrix' diagonal items.
     * @return a_i = A_i_i {1..n}
     */
    fun diagonal(): List<Double> = (0..Math.min(height, width) - 1).map { this[it, it] }

    /**
     * Calculates square matrix determinant.
     * @return det(this)
     */
    fun determinant(): Double {
        if (height != width)
            throw NotImplementedError()
        if (height == 1)
            return this[0, 0]
        if (height == 2)
            return this[0, 0] * this[1, 1] - this[1, 0] * this[0, 1]
        var result = 0.0
        for (i in rowIndices)
            result += (if (i % 2 == 0) 1.0 else -1.0) * this[i, 0] * minor(i, 0).determinant()
        return result
    }

    /**
     * Calculates the reversed matrix A^-1, so that A*A^-1 = E.
     * @return null, if this has determinant approximately equal to zero. this ^ -1 otherwise.
     */
    fun reversed(): Matrix? {
        val det = determinant()
        if (Math.abs(det) < REVERSE_DETERMINANT_EPSILON)
            return null
        return Matrix(
                rowIndices.map { i ->
                    colIndices.map { j ->
                        (if ((i + j) % 2 == 0) 1 else -1) * minor(i, j).determinant() / det
                    }
                })
    }

    /**
     * @return (A^T * A)^-1 * A^T
     */
    fun pseudoReversed() = transposed().let { (it * this).reversed()?.times(it) }

    /**
     * Makes a transposed matrix.
     * @return a_i_j = this_j_i
     */
    fun transposed(): Matrix = Matrix(
            colIndices.map { i ->
                rowIndices.map { j ->
                    this[j, i]
                }
            })

    /**
     * Makes a minor matrix, which is matrix without one row and one column.
     * @param a Row to delete
     * @param b Column to delete
     * @return Matrix of size [n-1 x n-1]
     */
    fun minor(a: Int, b: Int): Matrix = Matrix(
            (0..height - 2).map { i ->
                (0..width - 2).map { j ->
                    this[if (i < a) i else i + 1, if (j < b) j else j + 1]
                }
            })

    override fun toString() =
            StringBuilder().apply {
                for (i in rowIndices) {
                    for (j in colIndices) {
                        append(("%$TO_STRING_DIGITS.${TO_STRING_DIGITS}f").format(get(i, j)) + " ")
                    }
                    append("\n")
                }
            }.toString()

    /**
     * Returns a Vector-representation of specified row
     */
    fun getRow(i: Int): Vector {
        if (i < 0 || i >= rows.size)
            throw IllegalArgumentException("invalid row number\n")
        return Vector(ArrayList(rows[i]))
    }


    companion object {

        val REVERSE_DETERMINANT_EPSILON = 1e-14

        val SIMILARITY_EPSILON = 1e-6

        /**
         * Creates an unit matrix of size [n x n].
         * @param dimensions n -- size of the matrix.
         * @return A of size [n x n] with A_i_j = "i == j"
         */
        fun unit(dimensions: Int): Matrix = Matrix(
                (0..dimensions - 1).map { i ->
                    (0..dimensions - 1).map { j -> if (i == j) 1.0 else 0.0 }
                })


        fun zero(dimensions: Int): Matrix {
            return zero(dimensions, dimensions)
        }

        /**
         * Creates a matrix filled with zeros.
         * @param height Created matrix height.
         * @param width Created matrix width.
         * @return Matrix filled with zeros of size [height x width].
         */
        fun zero(height: Int, width: Int): Matrix = Matrix(
                (0..height - 1).map { i ->
                    (0..width - 1).map { j -> 0.0 }
                })

        val TO_STRING_DIGITS = 3

        fun byTrace(trace: DoubleArray): Matrix = Matrix(
                trace.indices.map { i ->
                    trace.indices.map { j ->
                        if (i == j) trace[i] else 0.0
                    }
                })

        fun parse(s: String): Matrix {
            val splitPoints = ArrayList<Int>()
            var bracesBalance = 0
            for (i in 0..s.length - 1) {
                if (s[i] == '{') {
                    ++bracesBalance
                } else if (s[i] == '}') {
                    --bracesBalance
                } else if (s[i] == ',' && bracesBalance == 0) {
                    splitPoints.add(i)
                }
            }
            splitPoints.add(s.length)
            val components = splitPoints.map { ArrayList<Double>() }
            for (i in splitPoints.indices) {
                val line = s.substring(if (i == 0) 0 else splitPoints[i - 1] + 1, splitPoints[i] - 1).trim { it <= ' ' }
                        .replace("[\\}\\{,]".toRegex(), " ").split("\\s+".toRegex()).dropLastWhile { it.isEmpty() }
                for (j in 1..line.size - 1)
                    components[i] += line[j].toDouble()
            }
            return Matrix(components)
        }
    }
}